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ABSTRACT 

Observations of very high speeds among pulsars in the Galactic disk present a puzzle regarding neutron 
stars in globular clusters. The inferred characteristic speed of single pulsars in the Galaxy is ~ 5 — 10 
times as large as the central escape speed from the most massive globular clusters. It is then reasonable 
Ch I to ask why any pulsars are seen in globular clusters, whereas, in fact, quite a large number have been 

^ ■ detected and as many as ^ 1000 are thought to be present in some of the richest clusters. A cluster that 

initially contains 10^ stars should be able to produce < 5000 NSs, simply based upon the mass function. 
, Therefore, it would seem that at least 10 — 20% of the NSs initially formed in a massive cluster must be 

retained until the current epoch. 

The recently derived distributions in natal "kick" speeds based upon isolated pulsars in the Galaxy 
are incompatible with the numbers of pulsars seen in certain globular clusters, if the cluster pulsars had 
isolated progenitor stars. This retention problem is a long-standing mystery. It has been suggested that 
the retention problem may be solved if one assumes that a large fraction of NSs in clusters were formed 
\^ • in binary systems. We present a thorough investigation of this possibility that involves a population 

' study of the formation and evolution of massive binary systems. 

We use a Monte Carlo approach to generate an ensemble of massive primordial binaries. Binary 
' component masses and orbital parameters are chosen at random from appropriate distribution functions. 

I Mass transfer begins when the more massive star evolves to fill its critical potential surface (Roche lobe). 

The mass transfer may be stable or dynamically unstable, depending on the structure of the mass donor 
I ' and the mass ratio of the components. Dynamically unstable mass transfer leads to a common-envelope 

O I phase and a dramatic reduction in the orbital separation, while a modest change in separation is expected 

■ if the transfer is stable. In either case, the orbital evolution is followed with an analytic prescription. It 

^ ' is assumed that the entire hydrogen-rich envelope of the initially more massive star is removed in the 

• • , mass transfer episode, exposing the star's helium core. The eventual collapse and supernova explosion 

^ ' of the core is accompanied by sudden mass loss and an impulsive kick to the newly- formed neutron star. 

k> I If we apply the large mean neutron star kick speeds inferred from pulsar observations, we find that most 

. binaries are unbound following the supernova, and all but a very small fraction of the liberated neutron 

' stars are ejected from the cluster. As expected, the majority of retained NSs have massive companions. 

These massive retained binaries are mostly the product of stable mass transfer, where the initially less 
massive star accretes a significant fraction of the envelope of the neutron star progenitor. Systems that 
undergo dynamically unstable mass transfer shrink dramatically and acquire large relative orbital speeds 
(typically > 200 km s~^). Sudden mass loss in the subsequent supernova explosion transforms > 10% 
of the orbital speed into translational motion. If, in addition, the newly-formed neutron star receives a 
large kick, it is likely that the system will escape from the cluster. 

Our "standard model" involving the formation of neutron stars in binary systems predicts that ~ 5% 
of the neutron stars initially formed in a massive cluster can be retained. Over a wide range of model 
parameters, the retention fraction varies from ^ 1 — 8%. When a number of other effects are taken into 
account, e.g., a reasonable binary fraction among massive stars, the retention fraction becomes several 
times smaller. Therefore, we suggest that perhaps the conventional thinking regarding neutron star kicks 
must be modified or that a new paradigm must be adopted for the evolution of some of the most massive 
globular clusters in the Galaxy. 

Subject headings: globular clusters: general — stars: neutron 



1. INTRODUCTION 

A growing body of observational and theoretical evi- 
dence suggests that some massive globular clusters may 



contain more than ^ 1000 neutron stars (NSs). However, 
the presence of even as few as ^ 100 NSs is difficult to 
reconcile with the large NS "kicks" inferred from proper 
motion studies of single, young radio pulsars in the Galac- 
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tic disk. The problem is that globular clusters have central 
escape speeds < 50 km s""'^, while it is widely thought that 
most NSs are born with speeds > 200 km s~^. This is the 
essence of the NS retention problem in globular clusters. 

If one accepts the conventional wisdom regarding NS 
kicks, then only a very small fraction of NSs that are 
remnants of isolated progenitors should be retained in a 
globular cluster. Hansen & Phinney (1997) found that a 
Maxwellian distribution in kick speeds, with a mean of 
^ 300 km s^^, is consistent with data on pulsar proper 
motions. This distribution predicts that only ~ 0.4% of 
NSs are born with speeds < 50 km s~^, and ~ 3% with 
speeds < 100 km s~^. If one adopts an initial mass func- 
tion derived from stars in the solar neighborhood (e.g., 
Kroupa, Tout. & Gilmore 1993), it can be shown that 

< 5000 NSs will be formed in a cluster that initially con- 
tains 10^ stars. A retention probability of 1% predicts that 

< 50 NSs should be present in a massive globular cluster 
such as 47 Tuc, where we have assumed that the cluster 
has not lost a significant fraction of its mass. Such a small 
number of NSs in 47 Tuc is incompatible with the obser- 
vational sample of more than 20 millisecond radio pulsars 
(Camilo et al. 2000) when selection effects are taken into 
account. 

Drukier (1996; using the results of Brandt & Podsiad- 
lowski 1995) and Davies & Hansen (1998) have demon- 
strated quantitatively that if a NS is formed in a massive 
binary systcnn, then there is a significant probability that 
the NS will remain in the binary following the supernova 
(SN) explosion, and that the recoil speed of the system 
could be sufficiently small to allow it to be retained in the 
cluster. While these studies provided a useful verification 
of the potential importance of massive binaries, they did 
not involve a systematic population study to determine a 
realistic NS retention fraction. 

Our primary goal in this paper is to make a detailed 
quantitative assessment of the role of massive binaries in 
retaining NSs in globular clusters. This calls for a real- 
istic description of the population of primordial binaries, 
as well as a sufficiently detailed consideration of the rele- 
vant stellar evolution processes that precede the first SN 
explosion. To this end, we have developed a Monte Carlo 
population synthesis code that follows each of an ensem- 
ble of massive, primordial binaries from the main-sequence 
phase, through any important episodes of mass transfer, 
up to and immediately beyond the time of the first SN. 

We have attempted to make the paper self-contained, 
so that much of the relevant background material and as- 
sociated references are provided. The paper is organized 
as follows. In § 2 we review the evidence indicating that 
NSs may be quite abundant in certain massive globular 
clusters. An historical overview of the debate regarding 
NS kicks is presented in § 3. In § 4 we describe the various 
elements of our Monte Carlo population synthesis code. A 
semi-analytic treatment of massive binary population syn- 
thesis and NS retention is given in § 5, which we intend to 
facilitate the interpretation of the results of our detailed 
numerical calculations presented in § 6. The main results 
of our study are reviewed in § 7, and we evaluate the pos- 
sibility that binaries provide a robust solution to the re- 
tention problem. Finally, we speculate in § 8 on possible 
alternative solutions to the NS retention problem. 



2. NEUTRON STARS IN GLOBULAR CLUSTERS 

Several tens of millisecond pulsars (MSPs), a dozen 
bright X-ray sources, and numerous low- luminosity X-ray 
sources have been detected in the Galactic globular cluster 
system. See Table 1 for a list of clusters that may contain 
large numbers of NSs. The nature of the pulsars is clear: 
these are rapidly spinning NSs, many of which have bi- 
nary companions. The luminous cluster X-ray sources are 
all low-mass X-ray binaries (LMXBs) powered by accre- 
tion onto a NS. An accepted familial relationship exists 
between LMXBs and the majority of MSPs, the former 
being the evolutionary progenitors of the latter. Recent 
observations (Grindlay et al. 2001) provide tantalizing ev- 
idence that many of the low- luminosity X-ray sources may 
be MSPs for which radio pulsations have not yet been de- 
tected (sec, however, Pfahl & Rappaport 2001). 

More refined pulsar searches, deeper X-ray observations, 
and thorough theoretical population studies will advance 
our understanding of the cluster population of NSs, and in 
turn may provide powerful new insights into the formation 
and evolution of globular clusters. We now briefly review 
what is known and what is speculated regarding NSs in 
globular clusters. 

2.1. Millisecond Pulsars 

Camilo et al. (2000; see also Freire et al. 2000), using 
the Parkes radio telescope, have detected more than 10 
MSPs in the globular cluster 47 Tuc, bringing the current 
total to over 20. With a cursory analysis of the selection ef- 
fects and a reasonable pulsar luminosity function, Camilo 
et al. (2000) estimated that 47 Tuc may contain 200 po- 
tentially observable MSPs, and therefore the total number 
of NSs in 47 Tuc is expected to be even larger. 

47 Tuc is a rich and relatively nearby cluster (distance 
~ 4.5 kpc), and thus has been the subject of much study, 
especially in recent years. Unfortunately, other clusters 
with properties similar to 47 Tuc have not yet received 
as much attention. However, there is compelling observa- 
tional evidence that at least one other cluster contains in 
excess of 100 NSs. Pruchter & Goss (2000) found signif- 
icant diffuse radio emission from the core of the massive 
globular cluster Terzan 5. These authors estimated that 
this diffuse component may be attributable to between 60 
and 200 potentially observable pulsars, and claim that this 
range may represent a rather severe underestimate of the 
total number of pulsars in the cluster. Thus far, two MSPs 
have been detected in Terzan 5 (Lyne, Mankelow, Bell, & 
Manchester 2000). 

2.2. X-ray Sources 

Many of the members of the well-known class of bright 
X-ray sources in globular clusters, with X-ray luminosities 
Lx ~ 10^^ - 10^^ ergs s~^ (see Deutsch, Margon, & An- 
derson 2000, and references therein), exhibit type I X-ray 
bursts (Lewin, van Paradijs, & Taam 1993) and are there- 
fore accreting NSs in binary systems. In fact, 7 of the 12 
known bright sources have a well-measured or constrained 
binary period (Deutsch, Margon, & Anderson 2000). Each 
of these objects resides in a different globular cluster (Ta- 
ble 1). While this sample of X-ray binaries does not con- 
stitute a large number of NSs in itself, the existence and 
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properties of these systems may have important impUca- 

tions regarding the cvohition of the NS population in their 
respective host clusters (see § 2.3 below). 

Not as well known, and certainly not as well under- 
stood, is the class of low-luminosity cluster X-ray sources 
(e.g., Johnston & Verbunt 1996; Verbunt 2001), with 
Lx ^ W"^^ — 10"^^ ergs s""'^, where the lower limit is set by 
detection sensitivities. Prior to the launch of the Chandra 
X-ray Observatory, fewer than 50 of these faint sources had 
been discovered in the entire Galactic globular cluster sys- 
tem (see Verbunt 2001 for a recent analysis of the ROSAT 
database), primarily with the ROSAT and Einstein satel- 
lites. However, recent deep Chandra observations of 47 
Tuc have revealed > 100 faint sources in this cluster alone 
(Grindlay et al. 2001), whereas only 9 had been confirmed 
previously (Verbunt & Hasinger 1998). There is growing 
evidence that the majority of the faint X-ray sources in 47 
Tuc may be NSs, perhaps MSPs that have not yet been 
detected at radio wavelengths. All of the 15 MSPs in 47 
Tuc with well-measured radio timing positions (Freire et 
al. 2000) have counterparts in the Chandra images (Grind- 
lay et al. 2001). Further multiwavelength observations are 
required to determine the true distribution of objects con- 
tributing to the population of low-luminosity sources. 



2.3. Theoretical Considerations 

Variations in the number of detected radio pulsars and 
X-ray sources from cluster to cluster may be attributed 
to the distances of the clusters, selection effects inher- 
ent in the observations, as well as differences between the 
intrinsic NS populations. Predictions of the total num- 
ber of potentially observable NSs - in the form of MSPs 
or accretion-powered X-ray sources - present in a glob- 
ular cluster are difficult. Empirical likelihood estimates 
based on the observational sample are hindered by small- 



number statistics and uncertainties regarding selection ef- 
fects. Theoretical studies aimed at accounting for the 
numbers and properties of the detected pulsars involve 
models that utilize various uncertain stellar evolution and 
dynamical processes. 

Large-scale population studies of the formation and evo- 
lution of X-ray binaries and MSPs in globular clusters have 
only recently been undertaken (see Davies 1995; Sigurds- 
son & Phinney 1995; Rasio, Pfahl, & Rappaport 2000; 
Rappaport et al. 2000). The dense stellar environment in 
a globular cluster allows for dynamical binary formation 
channels not possible in the Galactic disk, such as two- 
body tidal capture (e.g., Fabian, Pringle, & Rees 1975; 
Rasio & Shapiro 1991; DiStefano & Rappaport 1992), 
and three- and four-body exchange processes (e.g.. Hills 
1976; Hut, Murphy, & Verbunt 1991; Sigurdsson & Phin- 
ney 1993; Bacon, Sigurdsson, & Davies 1996; Rasio, Pfahl, 
& Rappaport 2000). The absolute probabilities of dy- 
namical encounters depend on the local stellar environ- 
ment and thus implicitly on the dynamical evolution of 
the cluster. This nonlinear linkage between local dynami- 
cal processes and the global cluster evolution poses signif- 
icant computational problems, but the potential rewards 
are far-reaching. Such population studies promise to be 
a powerful tool that relates the current NS population to 
the formation and evolution of globular clusters. 

Preliminary calculations (Rasio, Pfahl, & Rappaport 
2000; Rappaport et al. 2000) indicate that a large initial 
pool of single NSs (~ lO'') may be required to explain 
the handful of very short-period binary radio pulsars in 
47 Tuc; i.e., the formation eSiciency is quite low. Short- 
period binary MSPs and LMXBs in other clusters may be 
good indicators of an initially large number of NSs in those 
clusters as well. The purpose of the present paper is to in- 
vestigate the conditions that favor the retention of such a 
large NS population. 
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3. NEUTRON STAR KICKS 

It was suggested four decades ago by Blaauw (1961; see 
also Boersma 1961), before the discovery of the first radio 
pulsar (Hewish et al. 1968), that there may be a pop- 
ulation of massive-star remnants (neutron stars) in the 
Galaxy with anomalously large space velocities, acquired 
when the NS progenitor explodes in a binary system. If the 
orbit is circular and the explosion is spherically symmetric, 
then the compact remnant is liberated from its companion 
if more than one-half of the initial mass of the binary is 
lost in the SN. Blaauw (1961) proposed this mechanism to 
explain the famous "runaway" O and B stars (see Hooger- 
werf, de Bruijne, & dc Zceuw 2001 for a recent discussion). 

Gunn & Ostriker (1970) analyzed the data on some forty 
pulsars known at that time. They found that the pulsars 
had an average scale height of ^ 120 pc above the Galac- 
tic plane, in contrast to the smaller scale height (~ 50 pc) 
of massive O and B stars. Gunn & Ostriker (1970) con- 
cluded that pulsars are born with an average speed of 
~ 100 km s~^, and attributed this speed to a binary origin. 

By 1980, 26 pulsars had reasonably accurate interfero- 
metrically determined proper motions. Lyne, Anderson, 
& Salter (1982) showed that the observed pulsars had a 
three-dimensional rms speed of ^ 210 km s"-'^, and a ver- 
tical scale height of ~ 350 pc. Thus, pulsars are even faster 
and more dispersed than indicated in the study by Gunn 
& Ostriker (1970). Based on the data of Lyne, Ander- 
son, & Salter (1982), Paczyhski (1990) adopted a modified 
Lorentzian distribution of pulsar speeds. 
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where vq is the rms speed. The best agreement with the 
data was found with vq = 270 km s~^. 

The next major step in the study of the kinematics of 
the Galactic pulsar population came with the work of Lyne 
& Lorimer (1994). This analysis considered a total of 
99 pulsars, 86 of which had interferometric proper mo- 
tions or upper limits, and an additional 13 pulsars with 
only scintillation speed measurements (e.g., Cordes 1986). 
Helfand & Tadcmaru (1977) and Cordes (1986) noted a 
strong selection effect against detecting high-velocity pul- 
sars. Fast pulsars born in the Galactic disk rapidly move 
away from the plane and may only spend a small fraction 
of their radio lifetimes in the flux-limited volume. The 
vertical egress of fast pulsars thus implies a mean speed 
for the detectable pulsar population that is lower than the 
true mean pulsar birth speed. Lyne & Lorimer (1994) at- 
tempted to account for this effect by utilizing only those 
pulsars with characteristic ages less than 3 Myr, thus re- 
ducing the analyzed sample to 29 objects. Using a revised 
dispersion-measure distance scale (Taylor & Cordes 1993), 
Lyne & Lorimer (1994) derived a mean pulsar speed of 
^ 450 km s^^ and indicated that < 1% of pulsars have 
speeds less than 50 km s~^. 

The larger proportion of fast pulsars (speeds > 
200 km s""'^) began to cast doubt on the notion 
that the high velocities were acquired as a re- 
sult of the binary "slingshot" effect (see, however, 
Iben & Tutukov 1996 and the rebuttal by van den 
Heuvel & van Paradijs 1997). Shklovskii (1970) 
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Fig. 1. — The Maxwellian and modified Lorentzian kick distri- 
butions, sliown for iT(km s~^) = {50,100,200} and j;o(km s~^) = 
{100, 300, 500}, respectively. 

was the first to suggest that large pulsar speeds may result 

from an asymmetry intrinsic to the SN explosion. This is 
now a widely held view. Proposed physical mechanisms 
for natal NS kicks include purely hydrodynamical models, 
as well as scenarios that involve asymmetric neutrino emis- 
sion, possibly correlated with the orientation and strength 
of the magnetic field within the proto-ncutron star (Ar- 
ras & Lai 1999). For a recent review of various natal kick 
mechanisms, sec Lai (2000). Harrison & Tademaru (1975; 
see also Lai, Chernoff, & Cordes 2001) have proposed an 
alternative mechanism for the acceleration of pulsars that 
involves asymmetric dipole radiation after the pulsars has 
been formed, the so-called "electromagnetic rocket" mech- 
anism. 

Since the work of Lyne & Lorimer (1994), a number of 
authors have reanalyzed the pulsar proper motion data 
with more sophisticated treatments of selection effects 
(e.g., Lorimer, Bailes, & Harrison 1997; Hansen & Phin- 
ney 1997; Cordes & Chernoff 1998). Each of these stud- 
ies employed a different statistical methodology, but all 
found a mean three-dimensional pulsar speed in excess 
of 300 km s^^. Lorimer, Bailes, & Harrison (1997) and 
Hansen & Phinney (1997) found that a Maxwellian distri- 
bution in natal kick speeds, of the form 



P(i'fe) = 



IT 



(2) 



was reasonably consistent with the proper motion data, 
with a ^ 300 km s~^ and ^ 200 km s~^, respectively, 
in the two studies. Cordes & Chernoff (1998) found ev- 



idence for a bimodal kick distribution, with > 



of 



the pulsars contained in a Maxwellian component with a 
dispersion of cr < 200 km s~^, and the remaining NSs 

contained in a high-speed component with a dispersion of 



700 km 



In all of these investigations, a Maxwellian 



distribution was adopted by convention. However, both 
Lorimer, Bailes, & Harrison (1997) and Hansen & Phin- 
ney (1997) stress that, while it is likely that the mean birth 
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speed of pulsars is large, and that high- velocity pulsars are 
under-represented in the observational sample, the data do 
not constrain the functional form of the underlying kick 
distribution (see also. Fryer, Burrows, & Benz 1998), nor 
do the data suggest a physical mechanism for the kick. 

For completeness, we also list a number of pieces of in- 
direct evidence for NS kicks. 

(1) Pulsar velocities > 1000 km have been inferred 
from pulsar-supernova remnant associations (e.g., Caraveo 
1993; Frail, Goss, & Whitcoak 1994) by dividing the dis- 
tance from the rough geometrical center of the remnant by 
the characteristic spin-down age of the pulsar; however, 
the reliability of some of these estimates is questionable 
(see Gaensler & Frail 2000). 

(2) For two compact double neutron star systems not in 
a globular cluster, PSR 1913-hl6 and PSR 1534-hl2, Fryer 
& Kalogera (1997; see also Flannery & van den Heuvel 
1975) estimate that the most recently formed NS in each 
case had a kick of > 200 km s~"^. 

(3) Kaspi ct al. (1996) have found that the B-star com- 
panion to the pulsar PSR J0045-7319 in the SMC has a 
spin that is inclined with respect to the orbital angular 
momentum vector, and may, in fact, have retrograde rota- 
tion. This is convincing evidence that the pulsar received 
a kick with a significant component perpendicular to the 
pre-SN orbital plane. 

(4) Systems with moderate-to-high eccentricities (e > 
0.3) seem to dominate the sample of Bc/X-ray bina- 
ries with measured orbital parameters (see Bildsten et 
al. 1997). These eccentricities may require kick speeds 
> 50 km (Verbunt & van den Heuvel 1995). How- 
ever, a number of Bc/X-ray binaries (XTE J1543-568, 4U 
0352-F30/X Per, 2S 1553-54, GS 0834-43, 7 Gas; see Ta- 
ble 6) exhibit relatively low eccentricities (e < 0.2) and 
orbital periods sufficiently long that tidal circularization 
could not have played a role (Delgado-Marti et al. 2001). 
A detailed study of X Per/4U 0352+30 (Dclgado-Martf et 
al. 2001) revealed that the present orbit (-Porb = 250 d, 
e = 0.11) is entirely consistent with the neutron star hav- 
ing been born with no kick. This may indicate that the 
underlying NS kick distribution is not simply related to 
the speed distribution of isolated pulsars (see § 8.2). 

Large uncertainties regarding natal NS kicks imply that 
a comprehensive population study must incorporate a va- 
riety of kick distributions. In our simulations of the for- 
mation of NSs in globular clusters, we c;onsidcr a single- 
component Maxwellian kick distribution as well as the 
modified Lorentzian proposed by Paczynski (1990) (see 
Figure 1). 

4. POPULATION SYNTHESIS OF MASSIVE BINARIES 

Aside from what little information can be inferred from 
the current cluster population of stars, there is insuffi- 
cient data to tightly constrain the distributions in main- 
sequence masses and orbital parameters of the massive pri- 
mordial binaries in globular clusters. However, a combi- 
nation of theoretical and empirical evidence suggests that 
the process of star formation (single and binary) is quite 
generic in a qualitative sense (e.g., Elmegreen 2000, and 
references therein), and that the general "rules" that ap- 
ply to star formation in the Galactic disk may also apply 
to regions of high stellar density, such as globular clus- 



ters. Fortunately, there is a large body of work on the 
binary population in the Galactic disk (e.g., Abt & Levy 
1978; Kraicheva et al. 1978; Duquennoy & Mayor 1991). 
For these reasons, and for lack of any concrete alternative 
siiggestions, our population synthesis study of massive bi- 
naries in globular clusters parallels similar studies of mas- 
sive binaries in the disk (e.g., Podsiadlowski, Joss, & Hsu 
1992; Terman, Taam, & Savage 1998). 

Our population synthesis code is comprised of three ba- 
sic elements, corresponding to the three main evolution- 
ary stages of a massive binary: (i) the masses and orbital 
parameters are chosen from appropriate distribution func- 
tions, (ii) analytic approximations are used to follow the 
binary stellar evolution through any important episodes of 
mass transfer, and (iii) the dynamical influence of the first 
SN explosion is computed. We now discuss each of these 
elements in detail. 

4.1. Primordial Binaries 

We define a massive primordial binary as a system where 
at least one component has a large enough main-sequence 
mass that it would yield a neutron star remnant if left 
to evolve in isolation. The main-sequence mass threshold 
for NS formation is ~ 8 Mq, with a weak dependence on 
metallicity (see § 5). Hereafter, we refer to the initially 
more massive component of the binary as the primary and 
the initially less massive component as the secondary. The 
main-sequence masses of the primary and secondary are 
denoted by Mi and M2, respectively, and we define the 
main-sequence mass ratio as q = M2/M1 < 1. 

Various authors (e.g., Kraicheva ct al. 1978; Duquen- 
noy & Mayor 1991) have found that the primary masses 
in close binaries are consistent with the initial mass func- 
tion (IMF) of single stars. We draw primary masses from 
a single power-law IMF, which is appropriate for massive 
stars (sec Miller & Scalo 1979; Scalo 1986; Kroupa, Tout, 
& Gilmore 1993): 

p{M,) = {x- 1){M^^^ M,-^t^)-'Mr , (3) 
where A'/i ^in ^ 8 Mq. We choose values of x in the 
range 2-3, where x = 2.35 corresponds to a Salpeter IMF 
(Salpeter 1955). For our standard model (see § 6.2.1), we 
adopt X = 2.5. An isolated star with Mi > Mimax is 
assumed to leave a black hole remnant rather than a NS. 
The precise value of Mi^^ax is unknown, but is probably 
> 30 Mq. Because the IMF sharply decreases with in- 
creasing mass, our results do not depend strongly on the 

value of Mi^max- 

It is generally believed that the main-sequence masses 
of binary components are correlated (e.g., Abt & Levy 
1978; Garmany, Conti, & Massey 1980; Eggleton, Fitch- 
ett, & Tout 1989). Observations of massive binaries sug- 
gest that equal masses may be favored (e.g., Abt & Levy 
1978; Garmany, Conti, & Massey 1980), presumably as a 
result of the formation process. However, serious selection 
effects hamper the determination of the true mass ratio 
distribution, and other authors have found that low-mass 
companions to massive primaries appear to be more likely 
(e.g.. Mason et al. 1998; Preibisch et al. 2000). For the dis- 
tribution function, p{q), we consider both increasing and 
decreasing functions of q, as encapsulated by the power- 
law form 

p{q) = {l + y)q^ , (4) 
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for y > ~1. We take y = 0.0 for our standard model (see 
§ 6.2.1). 

It is often assumed in population studies of the sort we 
have undertaken that the orbits of primordial binaries are 
circular. However, since the process of binary formation is 
only poorly understood, there is no a priori justification for 
suggesting that massive, primordial binaries should have 
small eccentricities. The situation is less clear in globu- 
lar clusters, where dynamical interactions may influence 
binary formation (e.g., Price & Podsiadlowski 1995; Bon- 
nell, Bate, & Zinnecker 1998). However, it is expected that 
the binary will eventually circularize if the orbit is suffi- 
ciently compact that mass transfer will take place. From 
here on we only consider circular primordial binaries. 

There exists at least one important caveat to the cir- 
cularization assumption stated above. Eccentricities of 
~ 0.5 are seen among the wide, interacting W Cephei 
binaries, which consist of a massive red supcrgiant and an 
early-type B star (e.g., Cowley 1969; Cowley, Hutchings, 
& Popper 1977). The VV Cephei systems are the widest 
of the known massive, interacting binaries, with periods 
> 10 yr, and their significant eccentricities may indicate 
that not enough time has elapsed for the system to cir- 
cularize. The growth of the large red-supergiant envelope 
occurs over a relatively short timescale of < 10^ yr, while 
the circularization timescale varies as {a/Rf (e.g.. Hut 
1981), where a is the binary semimajor axis and R is the 
radius of the supcrgiant. Therefore, the components of 
the observed VV Cephei binaries have experienced strong 
tidal interactions only very recently, in terms of the evo- 
lutionary history of the more massive star. 

The binary separation, a, is chosen from a distribution 
that is uniform in the logarithm of a (Abt & Levy 1978; 
see, however, Duquennoy & Mayor 1991): 



p{a) = 



In 



(5) 



For given component masses, the lower limit, amin, is 
determined from the constraint that neither star over- 
flows its Roche lobe on the main sequence. The upper 
limit, ttmax, is somewhat arbitrary; in practice, we assume 
amax = 103 AU. 

4.2. Mass Transfer: Overview 

Owing to its larger mass, the primary will be the first 
star to leave the main sequence. The subsequent binary 
evolution depends on the size of the Roche lobe that 
surrounds the primary. For circular orbits, the volume- 
equivalent radius of the Roche lobe of the primary is well 
approximated by the formula due to Eggleton (1983): 



a 



0.49 



0.6-Fg2/3in(i_^g-i/3) 



(6) 



where a in this equation represents the constant orbital 
separation. The radius of the Roche lobe of the secondary, 
i?L2, is obtained by replacing q in eq. (6) with 1/q. 

Left to evolve in isolation, the primary would grow to 
a maximum radius of ^ 500 — 2000 Rq (the value de- 
pends sensitively on mass, metallicity, and especially as- 
sumptions about stellar winds). So, if the orbit is suffi- 
ciently compact (i?Li 20 AU), the primary may grow to 
fill its Roche lobe. The value of arLi is used to determine 



the evolutionary state of the primary when it begins to 
transfer matter through the inner Lagrange point. 

It is particularly important to distinguish between mass 
transfer that is dynamically stable (proceeding on the nu- 
clear or thermal timescale of the primary) and mass trans- 
fer that is dynamically unstable (proceeding on the dy- 
namical timescale of the primary). For the case where the 
mass donor is more massive than the accretor, dynamical 
instability is typically attributed to one of two root causes. 
If the star grows faster than its Roche lobe (or the Roche 
lobe shrinks faster than the star does), a phase of run- 
away mass transfer may ensue. In particular, stars with 
deep convective envelopes tend to expand in response to 
mass loss (e.g., Paczyriski & Sienkiewicz 1972; Hjellming 
& Webbink 1987), while the Roche lobe generally shrinks. 
Also, for systems with extreme mass ratios, the primary 
may not be able to achieve synchronous rotation with the 
orbit, causing the components to spiral together; this is 
the classic Darwin tidal instability (Darwin 1879; see also 
Hut 1981) . The evolutionary state of the primary when it 
fills its Roche lobe is a good indicator of the physical char- 
acter of the subsequent mass transfer and binary stellar 
evolution. 

Following Kippenhahn & Weigert (1966; see also Lauter- 
born 1970 and Podsiadlowski, Joss, & Hsu 1992), we dis- 
tinguish among three evolutionary phases of the primary 
at the onset of mass transfer. Case A evolution corre- 
sponds to core hydrogen-burning, Case B refers to the shell 
hydrogen-burning phase, but prior to central helium igni- 
tion, and case C evolution begins after helium is exhausted 
in the core. These three cases (sec Figs. 2 and 3) provide 
a rough framework for categorizing binary stellar evolu- 
tion during mass transfer. Of course, not all systems will 
undergo Roche lobe overfiow. A large fraction of binaries 
will be sufficiently wide that the primary and secondary 
evolve as isolated stars prior to the first SN. We refer to 
such detached configurations as case D. 

Of the three broad categories of mass transfer, case A 
evolution is potentially the most problematic. Perhaps 
the most likely outcome is a merger of the two stars fol- 
lowing a contact phase, wherein both stars fill their Roche 
lobes, leaving a massive single star (Pols 1994; Wellstein, 
Langer, & Braun 2001). For a recent detailed discussion 
of the complex evolution processes that occur during case 
A mass transfer and the possible outcomes, see Nelson & 
Eggleton (2001) and Wellstein, Langer, & Brami (2001). 
Fortunately, the subtleties of case A evolution may essen- 
tially be overlooked in the present study. The range in 
orbital separations admitted by case A is ^ 3 — 20 Rq, 
which comprises ^ 5 — 10% of the primordial binary pop- 
ulation. In addition, if the majority of case A systems 
merge as expected, then a detailed consideration of case A 
evolution is unnecessary, since our focus is on how binarity 
impacts the retention problem. In our population synthe- 
sis code, case A mass transfer is treated in precisely the 
same way as early case B mass transfer (see below). This 
treatment is certainly unrealistic, but it is a highly opti- 
mistic scenario in regard to the retention problem, and so 
functions to provide the max;imum retention fraction for 
the case A systems. 

Two important subcases comprise case B. Early case 
B (case Be) mass transfer occurs when the primary fills 
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Fig. 2. — Evolution of the radius as a function of time for a star 
of mass 10 Mq and metallicity Z = 0.001. The range of radii for 
each case of mass transfer is labeled. Note that the radius decreases 
by a factor of ~ 10 immediately following central helium-ignition. 



Fig. 3. — Evolution of the radius as a function of time for a 
star of mass 1-5 Mo and metallicity Z = 0.001. Helium ignites in 
the core while that star is evolving through the Hertzsprung Gap. 
Steady helium-burning proceeds rapidly, and the radial evolution 
is barely perturbed during this phase, hence the term "failed blue 
loop." This type of evolution permits early, or radiative, case C 
mass transfer. 



its Roche lobe as it evolves through the Hertzsprung Gap 
(subgiant branch). In this case, the envelope of the pri- 
mary is still mostly radiative and the mass transfer is 
thought to be initially stable for a wide range of mass 
ratios. Late case B (case Bi) mass transfer occurs when 
the primary fills its Roche lobe as it evolves up the first 
giant branch. Case B; mass transfer is characterized by 
a deep convective envelope, in which case it is likely that 
mass transfer will initially take place on the dynamical 
timescale of the primary, leading to a common-envelope 
(CE) phase (see § 4.3). The range in orbital separations of 
case B systems is ~ 20 — 1000 Rq, which contains 30 — 40% 
of the primordial binary population. 

For stars of mass < 12 — 15 Af©, core helium-burning 
is typically accompanied by a significant decrease in stel- 
lar radius (see Fig. 2), so that the primary cannot fill its 
Roche lobe during this time. After helium is exhausted in 
the core, the star develops a deep convective envelope and 
begins to ascend the Hayashi track. We refer to this phase 
as late case C (case C;). We assume that case C; mass 
transfer is dynamically unstable, as we do for case B; (see, 
however, Podsiadlowski, Joss, & Hsu 1992; Podsiadlowski 
et al. 1994). 

For stars more massive than ^ 12 — 15 Mq, stellar evolu- 
tion calculations do not provide a self-consistent physical 
picture that simultaneously accounts for observations of 
massive stars in high- and low-metallicity environments 
(e.g., the Milky Way and the SMC, respectively). One 
main problem is to explain the ratio of blue to red super- 
giants among massive stars, as a function of metallicity 
(see Langer & Maeder 1995). For these massive stars, it 
is possible for helium to ignite in the core while the star 



is traversing the Hertzsprung gap and is still mostly ra- 
diative. It is unclear how the radius behaves during the 
subsequent phase of core helium-burning. The radius may 
shrink somewhat, so that Roche lobe overflow is impossi- 
ble during this phase. However, the radius may continue 
to increase after a short plateau (a so-called "failed blue 
loop" in the HR diagram; see Fig. 3). Therefore, it is the- 
oretically possible for mass transfer to begin during core 
helium-burning, but the allowed range in stellar radii is 
sufficiently small that we neglect this possibility. Follow- 
ing core helium-exhaustion, the star still has a radiative 
envelope. Early case C (case Cg) mass transfer is possi- 
ble before the primary reaches the base of the asymptotic 
giant branch and begins to ascend the Hayashi track, at 
which point the star is mostly convective and falls into the 
case C; category. Cases Ce and C; account for 15% and 
~ 10% of the primordial binaries, respectively. 

Depending on the choice of amax, the fraction of bina- 
ries that remain detached prior to the first SN (case D) is 
between 30% and 60%. For solar metallicity, stars of mass 
> 20 — 25 Mq may shed their hydrogen-rich envelopes 
following core helium-exhaustion as a result of prodi- 
gious wind mass loss (Maeder 1992). However, for low- 
metallicity environments like globular clusters {Z < 0.001 
typically), stellar wind mass loss is probably only impor- 
tant above ~ 30 - 40 Mq (Maeder 1992), although it 
should be emphasized that such mass loss is quite uncer- 
tain, both theoretically and observationally. In our simu- 
lations, we assume that no mass is lost prior to the SN ex- 
plosion. Clearly, since case D binaries are weakly bound, 
they make a negligible contribution to the NS retention 
fraction when the kick speeds are large. 
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The mass of the primary and its radius at the onset of 
Roche lobe overflow are sufficient to determine the cate- 
gory of mass transfer. Using the fitting formulae of Hurley, 
Pols, & Tout (2000), we compute the radius of an isolated 
star with the mass of the primary at different stages of its 
evolution (e.g., main sequence, core helium-ignition, base 
of the asymptotic giant branch). The Roche lobe radius of 
the primary falls into some range and determines if mass 
transfer is categorized as case A, Be, Bi, Ce, C;, or D. 

It is worth adding a note here regarding the possible 
outcomes of contact evolution. Unless the initial mass ra- 
tio is close to unity, it is likely that the majority of case 
A, Be, and Cg systems will evolve into a contact config- 
uration (Pols 1994; Nelson & Eggleton 2001; WcUstein, 
Langor, & Braun 2001). Shortly after the primary first 
fills its Roche lobe, the mass transfer rate rises to values 
of order Mi/rth.i ~ 10"^ - 10"^ Mq yr'^ where Tth,i is 
the thermal timcscale of the primary. Most of the envelope 
of the primary is removed during this initial rapid phase. 
The secondary reacts to donated material on its own ther- 
mal timescale, and thus can only accrete in an equilibrium 
fashion if it is very nearly coeval with the primary at the 
onset of Roche lobe overflow (i.e., q ^ 1). If the secondary 
is essentially unevolved at the time the primary fills its 
Roche lobe, the transferred gas will fill up the Roche lobe 
of the secondary and thus initiate a contact phase. 

Some attempts have been made to follow the evolu- 
tion of both stars during the contac;t phase (sec Kahlcr 
1989 for an overview of some of the important issues), 
but since the problem requires certain three-dimensional 
hydrodynamical elements, no complete physical descrip- 
tion of this process has emerged. However, there are two 
distinct possible outcomes: (i) evolution during contact 
is dynamically stable, with some fraction of the material 
ejected from the system, or (ii) there is sufficient drag on 
the secondary owing to the extended common stellar 
envelope surrounding the system - that the binary com- 
ponents spiral together on a dynamical timescale, possi- 
bly resulting in the ejection of the common envelope or 
a merger of the two stars. Dynamical spiral-in cannot be 
the generic consequence of contact evolution, since then 
we would have great difficulty in explaining the observed 
long-period high-mass X-ray binaries in the Galaxy (see 
§ 8.2). 

We expect that there is some critical mass ratio, qait, 
that separates stable and dynamically unstable mass 
transfer. In our simulations, we typically choose ^crit = 0.5 
for all case Be and Ce systems, which implies that 50% of 
these binaries undergo dynamically unstable mass transfer 
for a flat distribution in initial mass ratios. 

4.3. Mass Transfer: Analytic Prescriptions 

Dynamically stable mass transfer can be reasonably well 
characterized by two dimensionless parameters, not neces- 
sarily flxed during the evolution. The secondary accretes a 
fraction /3 of the mass lost by the primary. Generally, /3 is a 
function of the component masses, the rate at which mass 
is transferred from the primary, and the evolutionary state 
of the secondary. The complementary mass fraction, 1 — /?, 
escapes the binary system, taking with it speciflc angular 
momentum a, in units of the orbital angular momentum 
per unit reduced mass, (GMfca)^/^, where Mb = Mi+ M2 



is the total binary mass. 

For a circular orbit, the orbital angular momentum is 
given by 



(7) 



Mo 



(10) 



Logarithmic differentiation yields (e.g., Rappaport, Joss, 
& Webbink 1982) 

SJ _ SMi 6M2 _ 1 (5M(, 1 5a 

where it is assumed that orbit remains circular during the 
evolution. From the deflnition of the capture fraction, /3, 
it is clear that 5M2 = -fid Mi and 5Mb = (1 - 13) 5 M-^. 
We neglect the coupling between the rotation of the Roche 
lobe-filling primary and the orbit. Therefore, variation in 
the orbital angular momentum is attributed solely to sys- 
temic mass loss, and it follows that 

5J = a{GMba)^/'^5Mb . (9) 
We consider one of two modes of angular momentum 
loss during stable mass transfer: mass that is lost from 
the system takes with it (1) a constant fraction of the spe- 
cific orbital angular momentum (constant a), or (2) the 
specific angular momentum of the secondary. Constant 
values of a and /3 lead to one possible analytic solution of 
eq. (8) (Podsiadlowski, Joss, & Hsu 1992) 

where 

Ci = 2q;(1 - /?) - 2 

C2^-2aQ-l) -2. (11) 

Primes on the masses and semimajor axis indicate the 
values after some amount of mass has been transferred. 
This solution assumes /? > 0. The solution for /? = 
(totally non-conservative mass transfer) is obtained by re- 
placing (M^/Afa)*^" with exp[2a(M{ - Mi)/M2]; clearly, 
M2 remains fixed in this case. Figure 4 illustrates the 
evolution of (o//a)i as a function of the fractional mass 
loss, AAfi/Mi = 1 — M[/Mi, from the primary. For our 
standard model, we use eq. 10 to evolve the orbit, with 
a = 1.5, a value characteristic of mass loss through the L2 
point, and B = 0.75. These values imply that {a'/a)i ~ 1 
for AMi/Mi ~ 0.8. 

A second analytic solution can be obtained for the case 
where /3 is constant and where matter lost from the system 
takes away the specific angular momentum of the accret- 
ing star (perhaps in the form of jets or an axisymmetric 
wind). In this scenario we have 

-1 / j^^' \ -2 / \ -2//3 

2 yMbJ 

Note that the orbit always expands in this case. Further- 
more, it can be shown that eq. (12) is a weak function of 
/3, and so (a'/a)2 closely follows the /? = 1 solution for a 
wide range of initial mass ratios (within a factor of two). 

Mass that is removed from the primary on a dynamical 
timescale cannot be assimilated by the secondary, which 
can only accept matter on its much longer thermal read- 
justment timcscale. Consequently, the transferred mate- 
rial fills the Roche lobe of the secondary and a common- 
envelope (CE) phase is initiated (Paczynski & Sienkiewicz 
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\M2 



(12) 
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Fig. 4. Curves representing the dependence of the final-to- 
initial orbital separation, (a'/a)i, on the fractional mass loss from 
the primary, AMi/Afi, shown for ten different values of the mass 
capture fraction /3. The initial mass ratio was chosen to be g = 2/3 
and the dimensionless angular momentum-loss parameter was set 
to a = 1.5, a value characteristic of mass loss through the L2 point. 
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common-envelope evolution as a function of the initial mass ratio, 
q. The function is plotted for four different values of the fractional 
envelope mass, Me/Mi, of the primary. For this plot, we have 
chosen rjcE = 1-0 and A = 0.5. 



1972). As a result of hydro dynamic drag, the secondary 
spirals in toward the core of the primary, depositing a 
fraction rycE ^ 1 of the initial orbital binding energy 
into the CE as frictional luminosity (e.g., Meyer & Meyer- 
Hofmeister 1979; Sandquist, Taam, & Burkert 2000). For 
the binding energy of the envelope, we use the recent cal- 
culations of Dewi & Tauris (2000). If sufficient energy is 
available to unbind the envelope, what remains is a com- 
pact binary consisting of the secondary, which we assume 
is unaltered during the spiral-in process, and the helium- 
rich core of the primary. On the other hand, if the CE 
remains gravitationally bound to the system, drag forces 
will perpetuate the spiral-in and the two stars will merge. 

A simple energy relation determines the outcome of 
the CE and spiral- in (e.g., Webbink 1984; Dewi & Tau- 
ris 2000): 



GMiMe 

Xarhi 



VGE 



GM^Ma _^ GM1M2 



2a' 



2a 



(13) 



where Mf, and Mc are the envelope-mass and the core-mass 

of the primary, respectively. The left-hand side of eq. (13) 
is the envelope binding energy, where A is the structure 
constant computed by Dewi & Tauris (2000); for massive 
stars A ~ 0.5. Solving eq. (13) for a'/a, we find 



M,M2 
JVh 



M2 + T, \ • (14) 

If the secondary fills its Roche lobe for the computed fi- 
nal orbital separation, a', this effectively indicates that 
insufficient energy was available to unbind the CE, and 
we assume that a merger is the result. Figure 5 illustrates 
the dependence of (a'/a)cE on the initial mass ratio, for 
four values of the initial fractional envelope mass, Mg/Mi. 



In all cases where a stellar merger is avoided, we assume 
that the entire hydrogen-rich envelope of the primary is 
removed, either transferred stably through the inner La- 
grange point or expelled during a CE phase. By the time 
the primary reaches the base of the first giant branch (be- 
ginning of case B; evolution) its core is well developed, 
with a mass given approximately by (Hurley, Pols, & Tout 
2000) 



Mc ^ 0.1 



1.35 



(15) 



We assume that this is the mass of the helium core imme- 
diately following case Bg mass transfer as well, although 
it is expected to be somewhat smaller, since mass trans- 
fer interrupts the evolution of the primary. For case C 
and D evolution, the mass of the core may be larger by 
~ 0.5 — 1 Mq as a result of shell nuclear-burning. 

Following case B mass transfer, the exposed core of the 
primary is a nascent helium star. Helium stars with mass 
> 5 Mq probably experience mass loss in the form of a 
Wolf-Rayet wind, where the timescale for mass loss is com- 
parable to the evolutionary timescale (~ 10^ yr) of the star 
(see Langer 1989). Stellar winds are less important for he- 
lium stars of lower mass. However, if M^ < 3 Mq the 
star may grow to giant dimensions (Habets 1986b) upon 
central helium-exhaustion and possibly fill its Roche lobe, 
initiating a phase of case BB mass transfer (De Greve & de 
Loore 1977; Delgado & Thomas 1981; Habets 1986a). For 
the maximum radius of a helium star, we adopt a slightly 
modified version of the fitting formula derived by Kalogera 
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& Webbink (1998): 

log(-RHc,max/-R©) = 

r 2.3 Me < 2.5 

\o.O57[log(Mc/M0) - 0.17]-2-5 Mc > 2.5 . ^ 

This relation is consistent with the resuhs of Habets 
(1986b). It is expected that < 1 Mq is transferred to 
the secondary during the case BB phase. As an exam- 
ple, suppose the secondary has a mass of 10 Mq following 
case Be mass transfer, and the mass of the helium star is 
3 Mq, then the conservative transfer (/? = 1) of 0.5 Mq 
from the helium star leads to a ~ 30% expansion of the 
orbit. We consider case BB mass transfer by assuming 
that a fixed amount of mass (e.g., 0.5 Mq) is transferred 
conservatively to the secondary, and we expand the orbit 
accordingly. The inclusion of this process does not signif- 
icantly influence our results. 

4.4. Supernova Explosion 

At the end of the mass transfer phase, the result may 
be a stellar merger or a binary consisting of the secondary 
and the core of the primary. Subsequently, the remaining 
nuclear fuel in the primary core is consumed, leading to 
core-collapse and a SN explosion. The post-SN orbital pa- 
rameters are computed by taking into account the mass 
lost from the primary and the kick delivered to the newly- 
formed NS. In our simulations, we neglect the effect of the 
SN blast wave on the secondary. Statistically speaking, 
this assumption is well-justified, since only a small frac- 
tion of the binaries that we consider are sufficiently com- 
pact for the blast wave interaction to be important (e.g.. 
Wheeler, Lecar, & McKee 1975; Fryxell & Arnctt 1981; 
Livne, Tuchman, & Wheeler 1992; Marietta, Burrows, & 
Fryxell 2000). 

The magnitude and direction of the kick to the NS are 
often assumed to be uncorrelated. There is, as yet, no clear 
observational indication that the directions of NS kicks are 
preferentially aligned with respect to the spin of the NS 
progenitor. We assume that the orientations of the kicks 
are distributed isotropically. If the directions of the kicks 
are confined to a small cone perpendicular to the pre-SN 
orbital plane, then the net retention fraction is actually 
likely to be smaller than in the isotropic case, when the 
characteristic kick speed is larger than the typical pre-SN 
binary orbital speed (see Appendix A). 

One of two distinct outcomes follows the explosion: (i) 
the NS and the secondary are gravitationally bound, with 
new orbital parameters and a new CM velocity, or (ii) the 
binary is disrupted, with the NS and secondary receding 
along hyperbolic trajectories relative to the new CM. In 
the former case, the CM speed of the binary is compared 
to the cluster escape speed in order to determine if the 
NS, along with its binary companion, is retained in the 
cluster, while in the latter case the speed of the NS at in- 
finity, computed in the pre-SN CM frame of reference, is 
compared to the escape speed. 

When the characteristic kick speed is large, the mean 
eccentricity for the bound post-SN binaries may approach 
~ 0.7. In a significant fraction of these systems, the clos- 
est approach of the NS immediately after the SN is smaller 
than the radius of the secondary, with the likely outcome 
that the NS spirals in to the envelope of the secondary to 



form a Thorne-Zytkow object (Leonard, Hills, & Dewey 
1994; Brandt & Podsiadlowski 1995; see also § 6.5 for a 
more thorough discussion). 

The computational procedure that we use to compute 
the dynamical influence of the SN explosion on the binary 
system is outlined in Appendix B. This approach allows for 
the possibility that the pre-SN binary is eccentric, and, for 
completeness, the mathematical formalism also includes 
the effects of the interaction between the secondary and 
the SN eject a. 

In the majority of our simulations, we apply a single es- 
cape speed to all of the stars and binaries in question. If 
this escape speed is identified with the depth of the poten- 
tial well at the cluster center, then the computed retention 
fraction will be a maximum for a given set of parameters 
that describes the formation and evolution of the ensemble 
of binaries prior to the first SN. A more realistic cluster 
potential and spatial distribution of stars can only result 
in a smaller retention fraction. The flexibility of our pop- 
ulation synthesis code makes it straightforward to assess 
how a more realistic model of the cluster effects our results 
(see § 6.4). 

5. ANALYTIC ESTIMATES 

In order to facilitate the interpretation of our detailed 
population synthesis calculations, it is useful to develop 
some quantitative insight regarding how different assump- 
tions influence the NS retention fraction. This section is 
devoted to a semi-analytic population study of massive bi- 
naries and NS formation. We highlight the most profitable 
channels for retaining NSs in globular clusters. 

The IMF sets an upper limit to the total number of NSs 
that could have been formed in a cluster with some as- 
sumed initial total number of stars. For single stars with 
solar metallicity, the lowest initial stellar mass that yields 
a NS remnant is thought to be ~ 8 Mq. Globular clus- 
ters have a metal content that is typically < 10% of the 
solar value, and it is expected that the mass threshold for 
NS formation may be as low as ~ 6 Mq, but is probably 
not less than ~ 5 Mq (e.g., Marigo, Girardi, Chiosi, & 
Wood 2001, and references therein). The IMF of Kroupa, 
Tout, & Gilmore (1993) therefore predicts that between 
0.2% and 0.5% of single cluster stars should be sufficiently 
massive to produce NSs. If it is assumed that all massive 
stars are single, then a cluster that initially contains 10^ 
stars should produce < 5000 NSs. The situation is less 
clear for NSs born in binary systems, owing to the effects 
of mass transfer. 

Mass transfer in the case A and case Bg scenarios has 
the effect of interrupting the growth of the core of the 
primary, which implies a smaller final core-mass than if 
the primary evolved in isolation (e.g., Wellstein & Langer 
1999). Consequently, the stellar mass threshold for NS 
formation is increased above the conventional single-star 
value of 8 Mq , at least for solar metallicity. If mass trans- 
fer begins at a later stage of evolution (i.e., after central 
helium- burning), the final core-mass has already been es- 
tablished (to within < 1 Mq) and the mass threshold is 
basically unchanged. Thus, the fraction of primordial bi- 
naries that produce NS remnants is probably somewhat 
less than the corresponding fraction of single stars, for the 
same metallicity. The actual number depends on the pro- 
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portion of case A and case Be systems and on the relation 
between the final core-mass and the evolutionary state of 
the primary at the time it overflows its Roche lobe. 

Cases A, B, C, and D comprise roughly 5%, 25%, 25%, 
and 45%, respectively, of the primordial binary popula- 
tion. The 25% of case B systems are divided into ~ 20% 
case Be and ~ 5% case B;. Likewise, case C is comprised 
of ~ 15% case Cg and ~ 10% case C;. Depending on vari- 
ous assumptions, any of these percentages may increase or 
decrease by at most roughly one-half of tlic^ values given 
above. It is expected that ~ 35% of the primordial bina- 
ries evolve according to the case Be or Ce scenario; this 
has important implications for the retention problem. 

Stable and quasi-conservative mass transfer, which is 
expected in a significant fraction (depending on the value 
of 9crit) of the case Be and case Ce systems, has two no- 
table consequences: (i) the secondary accretes much of 
the hydrogen-rich envelope of the primary, and (ii) the fi- 
nal orbital separation is typically within a factor of a few 
of the initial separation (see Fig. 4). The increased mass 
of the secondary provides a deeper gravitational potential 
well for the helium star, which raises the likelihood that 
the orbit will remain bound following the SN and hence 
that the NS will be retained in the cluster. The modest 
change in orbital separation during stable mass transfer is 
in sharp contrast to the dramatic shrinkage that accom- 
panies CE evolution (see Fig. 5). For the case of small 
natal NS kicks, there is a clear dynamical distinction be- 
tween stable and unstable mass transfer with regard to 
the subsequent SN explosion. As a general rule, a NS kick 
can be considered "small" if it is appreciably less than the 
relative orbital speed of the components (see Brandt & 
Podsiadlowski 1995). 

First, consider the case of circular pre-SN orbits and 
vanishing kicks. If the binary is intact after the SN, reten- 
tion in the cluster is determined by the new center-of-mass 
speed, v'q^ (e.g., Blaauw 1961; Dewey & Cordes 1987): 

AMi 



VCM = 



-vi , 



(17) 



Mb - AMi 

where AMi is the mass lost in explosion (envelope of the 
primary) and vi is the pre-SN orbital speed of the pri- 
mary (helium star) about the CM. The mass Mf, and the 
orbital speed vi used above correspond to the conditions 
immediately before the explosion, and will generally differ 
from their main-sequence values as a resiilt of mass trans- 
fer. The factor that multiplies vi can be identified as the 
post-SN orbital eccentricity, e', which gives the memorable 
result, VQyi = e'vi. Disruption of the binary must occur 
if the mass lost in the explosion is more than one-half of 
the initial systemic mass. In this case it can be shown that 
the speed of the liberated NS is simply equal to vi . There- 
fore, regardless of whether or not the binary is unbound 
following the SN, the relevant speed that determines if the 
NS is retained in a globular cluster is proportional to the 
pre-SN orbital speed of the primary. 

These arguments are made more quantitative by consid- 
ering a prototypical binary with a range of initial separa- 
tions. Suppose this model primordial binary consists of a 
10 Mq primary and a 6 Mq secondary. A typical initial 
orbital separation for case Bg systems is ~ 0.5 AU. Im- 
mediately after a phase of stable mass transfer, the binary 
will consist of the ~ 2.2 Mq helium core of the primary 



and the - 10 - 14 Mq secondary, for /3 > 0.6. The bi- 
nary separation at this point is likely to be in the range 
~ 0.2 — 1.5 AU (see Fig. 4), in which case the helium 
core has an orbital speed of < 200 km s~^. A spher- 
ically symmetric SN will leave the system bound, with 
e' ~ 0.06 and v'q^ < 15 km s~^. There is then a good 
chance that such a binary would be retained in a globular 
cluster. For the case of CE systems that avoid a merger 
(most case B/ and C/ binaries), a typical initial separation 
is 5 AU. Following the expidsion of the CE, the essen- 
tially unaltered secondary orbits the core of the primary 
with a separation of < 0.05 AU (assuming hundred-fold 
decrease; see Fig. 5), giving the helium core an orbital 
speed of ~ 280 km s~^. For a 6 Mq secondary, the binary 
remains bound after a symmetric SN, but acquires a large 
recoil speed of ~ 30 km s~^. Thus, even in the case of 
small kicks, it is expected that a significant fraction of the 
post-CE binaries would be ejected from a typical globular 
cluster. 
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Fig. 6. — Retention probability as a function of the pre-SN orbital 
separation using the semianalytic formalism outlined in Appendix 
A. Curves arc shown for secondary masses M2 = 12 Mq (solid) 
and M2 = 6 Mq (dashed) and five different masses, Mc(Mq) = 
{2, 3, 4, 5, 6}, of the helium star prior to the supernova explosion. 
We have assumed an escape speed of 50 km s^^ and Maxwellian 
kicks with a = 200 km s~^. 



We now extend this discussion and allow for a distri- 
bution in NS kick speeds. The semi-analytic formalism 
that we employ is described in Appendix A, and the main 
results are displayed in Figure 6. Figure 6 is a plot of 
the retained percentage of bound NS binaries as a func- 
tion of the orbital separation immediately prior to the SN, 
where the escape speed is taken to be Ucsc ~ 50 km 
and kick speeds are distributed as a Maxwellian with 
a = 200 km s~^. The results are displayed for two dif- 
ferent secondary masses, M2 = 6 and 12 Mq, which are 
representative values for systems that undergo dynami- 
cally unstable mass transfer and stable mass transfer, re- 
spectively. A range of helium star masses, from Mc = 2 to 
6 Mq, is also considered. 
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There are a number of notable features in Fig. 6. First, 
the retention fractions are clearly larger for M2 = 12 Mq 
and a given core mass, as expected. Also, for a given M2 
and Mc, there is a maximum retention fraction located 
at some value of the initial orbital separation. The fall- 
off at large separations results from the high characteristic 
kick speed relative to the comparatively low orbital speeds; 
many of these systems are left unbound following the SN. 
At sufficiently small orbital separations, which correspond 
to high orbital speeds, the mean CM speed of the bound 
post-SN binaries exceeds the cluster escape speed, and 
thus accounts for the decreased retention fraction at small 
separations. The location and height of the peak depend 
on the helium star mass. In light of the discussion above 
regarding vanishing kicks, it is clear that a more massive 
helium star will result in a larger dynamical perturbation 
to the system at the time of the SN, simply by virtue of 
the increased mass loss. More mass loss results in a larger 
fraction of the orbital speed being transformed into CM 
speed for a bound post-SN binary (sec eq. [17]), and also 
raises the likelihood that the system will be disrupted if 
the characteristic kick speed is large. The combination of 
these effects explains the trends in Fig. 6, namely that, 
for a larger helium star mass, the height of the peak is 
reduced and its location shifts to larger separations (i.e., 
smaller orbital speeds). 

If we know the typical pre-SN component masses and or- 
bital separations among the systems that undergo stable 
or dynamically unstable mass transfer. Fig. 6 can be used 
to estimate the fraction of NSs that both remain bound 
to their companions following the SN and are retained in 
the cluster. The secondary masses used for Fig. 6 have 
already been appropriately chosen for this purpose. A typ- 
ical helium core mass is likely to be ~ 3 Mq. For the 
characteristic orbital separations, we chose 0.5 AU for the 
stable systems and 0.05 AU for the unstable systems (sec 
discussion above) . Restricting ourselves to the case B and 
C binaries, we should expect the ratio of stable to unsta- 
ble systems to be roughly 2 : 3. Now, reading numbers 
directly from Fig. 6, we estimate the percentage NS bi- 
naries retained following case B or C mass transfer to be 
~ 7%, with ~ 6% being systems where the mass transfer 
was stable and the remaining ~ 1% corresponding to the 
unstable systems. Of course, since case B and C systems 
comprise only ^ 50% of the primordial binary population, 
the net NS retention fraction is ^ 3.5%. This estimate is in 
accord with the results of the "standard model" discussed 
in the next section, where we calculate a net retention frac- 
tion of ^ 5% from all binary channels (that is, unweighted 
by the binary fraction among stars in the cluster). 

6. NEUTRON STAR RETENTION FRACTION 

6.1. Computational Procedure 

We have written two versions of our population synthe- 
sis code, each with same basic engine. One version applies 
a given distribution in kick speeds (e.g., a Maxwellian) and 
a single central escape speed. With this version of the code 
we are able to discern which individual evolutionary path- 
ways contribute most to the net NS retention fraction, and 
we can assess in a very detailed manner the influence of 
varying certain parameters. Only < 10^ primordial bina- 
ries are sufficient to obtain reliable statistics for the ~ 70 



distinct evolutionary channels followed in our code. 

The second version of our code is more global; here we 
are interested only in the net retention fraction. A large 
regular grid of kick speeds and escape speeds is established. 
We consider a range in escape speeds from to 100 km s""'^ 
and a range in kick speeds from to 1000 km s~^. At 
each position in the grid, an ensemble of binaries (typi- 
cally 2 X 10^) is generated and evolved. The output of 
the calculation is an array of retention fractions. It is then 
trivial to convolve the results with any of the kick distribu- 
tions discussed in § 3. For each escape speed and a given 
kick distribution, we compute a net NS retention fraction. 
An example of these retention curves is shown in Figure 
13. Additionally, we may also convolve the grid with a dis- 
tribution of escape speeds in order to gauge the influence 
of a realistic cluster potential (§ 6.4). 

6.2. Results 

6.2.1. Maxwellian Kicks 

We begin the discussion of our detailed population syn- 
thesis calculations by considering the Maxwellian kick dis- 
tribution. ThroTighout the rest of the paper we will refer 
to the "standard model." This reference model utilizes the 
Maxwellian kick distribution with a = 200 km s~"^ and a 
central escape speed of Vesc = 50 km s~^. The parame- 
ters that describe the primordial binary population and 
mass transfer for the standard model are listed in Table 3, 
model 5. 

Distributions of the primordial binary parameters for 
systems that undergo case B or C mass transfer are shown 
in Figure 7. We have not included systems that merge fol- 
lowing mass transfer, and hence the distribution in logo is 
not flat over the range shown, as would be expected from 
eq. 5. Furthermore, case A binaries are not included due 
to the small number of systems as well as the large uncer- 
tainties regarding their evolution. Figure 8 illustrates the 
correlation between the orbital separation and secondary 
mass for the systems in Fig. 7. The binaries that undergo 
dynamically unstable mass transfer and avoid a merger 
have a low mean secondary mass (~ 5 Mq) and a large 
mean initial separation (~ 5 AU). 

Figure 9 shows histograms of the binary parameters fol- 
lowing mass transfer for precisely the same systems in Fig. 
7. The distribution of secondary masses (sum of stable and 
unstable systems) shows a clear bimodality, with peaks 
around 5 and 12 Mq. This is simply a consequence of 
the distinction between stable and dynamically unstable 
mass transfer. If the mass transfer is stable, a secondary 
of mass > 4 Mq (assuming ^crit ^ 2) accretes a substantial 
amount of material, while the secondary mass is assiimed 
to be unchanged if the system evolves through a CE phase. 
Therefore, the peak in the distribution of secondary masses 
for the stable systems is shifted to a higher value, and the 
distribution as a whole is broadened, thus reducing the 
height of the peak relative to the secondary mass distri- 
bution for the unstable systems. Also noteworthy is the 
broad peak in orbital separations (summed distribution) 
centered at ~ 0.1 AU, which results from the overlap of 
the stable and unstable systems. 

A scatter plot of the secondary mass and orbital separa- 
tion following mass transfer is shown in Figure 10. There 
is a clearly-deflned boundary that marks the Roche lobe 
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Fig. 7. — Distributions of masses and orbital separations of pri- 
mordial binaries that undergo case B or case C mass transfer and 
which do not merge. Hatched regions indicate systems that un- 
dergo stable mass transfer (+45°) and dynamically unstable mass 
transfer (-45°). The histogram that encloses the hatched regions 
is the sum of the distributions in stable and unstable systems. 



Fig. 8. — Scatter plot of circulajrized orbital separation versus 
secondary mass for the primordial binaries shown in Fig. 7. Filled 
and unfilled circles indicate systems that undergo stable and dy- 
namically unstable mass transfer, respectively. 
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Fig. 9. — Distributions of masses and orbital parameters of sys- 
tems that have undergone case B or case C mass transfer and which 
have not merged. The hatchings have the same meaning as in Fig. 
7. These parameters indicate the state immediately prior to the 
supernova explosion of the helium core of the primary. 
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Fig. 10. — Scatter plot of separation versus secondary mass for 
the post-mass transfer binaries shown in Fig. 9. Note the well- 
defined boundary at small separations that marks the Roche lobe 
of the secondary. Filled and unfilled circles indicate systems that 
have undergone stable and dynamically unstable mass transfer, re- 
spectively. 
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Fig. 11. — Distributions of binary parameters of systems that 
have undergone case B or C mass transfer, have been left bound 
following the supernova explosion, and have been retained in the 
cluster. Note that many of the systems with log(a/AU) < —1.3 will 
experience a coalescence of the NS and the secondary immediately 
following the SN. 



Fig. 12. — Scatter plot of pcriastron separation versus secondary 
mass for the same systems as in Fig. 11. The open and filled circles 
have the same meaning as in Figs. 8 and 10. 
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Fig. 13. — Percentage of NSs retained as a function of the cluster 
escape speed, for a MaxwcUian distribution of kick sjjccds. Heavy 
curves correspond to binary channels, while the lighter curves are 
for single stars only. 



Fig. 14. — Percentage of NSs retained as a function of the clus- 
ter escape speed, for a modified Lorcntzian distribution of kick 
speeds. Heavy curves correspond to binary channels, while the 
lighter curves are for single stars only. 
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Table 2 

Branching Percentages for Maxwellian Kicks 



Outcome A B, C/ D 



Total 5.40 21.51 3.11 15.03 6.69 48.27 

SlabloMT* 2.73 10.83 0.00 7.56 0.00 0.00 

Unstable MT (CE"^) 2.68 10.67 3.11 7.47 6.69 0.00 

Merger Following MT 3.19 10.64 2.08 5.82 0.30 0.00 



a = 200 km 

Unbound following SN 0.78 7.52 0.47 7.42 3.96 48.26 

Bound following SN 1.43 3.34 0.56 1.78 2.44 0.01 

Merger following SN 0.53 0.81 0.24 0.46 0.93 0.00 

Retained single NS 0.01 0.10 0.00 0.08 0.04 0.22 

Retained binary NS 0.93 2.60 0.19 0.96 0.67 0.01 

Total retained 0.94 2.70 0.19 1.03 0.71 0.23 



a = 100 km s^^ 

Unbound following SN 0.18 4.50 0.17 5.30 2.05 48.22 

Bound following SN 2.03 6.36 0.86 3.91 4.34 0.05 

Merger following SN 0..33 0.56 0.19 0.38 0.92 0.00 

Retained single NS 0.01 0.27 0.01 0.31 0.09 1.54 

Retained binary NS 1.52 5.71 0.63 2.90 2.46 0.05 

Total retained 1.53 5.98 0.64 3.21 2.55 1.60 



a = 50 km s ^ 

Unbound following SN 0.00 1.73 0.01 2.76 0.63 47.91 

Bound following SN 2.20 9.13 1.02 6.45 5.77 0.36 

Merger following SN 0.07 0.13 0.07 0.10 0.46 0.00 

Retained single NS 0.00 0.31 0.00 0.55 0.08 9.90 

Retained binary NS 1.69 8.49 0.94 5.35 4.41 0.36 

Total retained 1.69 8.80 0.94 5.90 4.49 10.26 



=MT = Mass Transfer. 
''CE = Common-Envelope. 



radius of the secondary for a given M2 and M^- Systems 
to the loft of this boundary have merged following mass 
transfer, where the majority of merged binaries result from 
dynamically unstable case Be and Ce mass transfer. From 
Table 2 we see that roughly one-half of the case Be and 
Ce systems are expected to merge following mass trans- 
fer. This factor of one-half is a direct consequence of our 
choice of gcrit = 1/2. Nearly 50% more stable systems re- 
sult if we set gcrit = 1/4, but additional dilution factors 
lead to a net retention fraction that is only ~ 30% larger 
- that is, ~ 1.3 times as large - when v^sc = 50 km s~^ 
and a = 200 km s^^ (see Table 3). 

Figures 11 and 12 show distributions relevant to the 
bound and retained binaries immediately following the SN 
that underwent case B or C mass transfer (a subset of the 
binaries in Figs. 7 and 9). Small pcriastron separations 
(< 1 AU) among the retained NS binaries indicate that the 
secondaries in most of these systems, the majority of which 
have a mass > 10 Mq , will transfer material to the NS at 
some stage; in fact, in some cases (log(a/AU) < —1.3) the 
radius of the secondary is larger than the pcriastron sepa- 
ration immediately after the SN, indicating an immediate 
coalescence. Each of these points (mass transfer and coa- 
lescence) is discussed in § 6.5. Also, note that the speed 
distribution of the retained binaries has significant values 
all the way up to the escape speed. A more realistic clus- 
ter potential and spatial distribution of stars is therefore 
likely to result in a marked decrease in the net retention 
fraction, since the fastest of the binaries in Fig. 11 would 
be preferentially removed removed from the retained pop- 



ulation (see § 6.4). 

Table 2 and Figure 13 are the main results of our reten- 
tion study for Maxwellian kicks. The importance of case 
Be and Ce systems, which contribute a large number of 
bound and retained binaries (see § 5), is clear in Table 
2. Figure 13 shows the percentage of NSs retained in a 
cluster as a function of the central escape speed (applied 
to all stars and binaries); the curves are not weighted by 
the binary fraction. For a — 200 km s~^ the retention 
fraction is as large as ~ 2% for single stars and ~ 10% 
for binaries when v^sc = 100 km s~^. It is evident from 
Table 2 and Fig. 13 that the retention problem is elimi- 
nated for cr = 50 km s~^, where the retained fraction of 
NSs with isolated progenitors is ^ 10%, roughly one- half 
of the binary contribution. 

Finally, in order to gauge how the net retention fraction 
changes when the free parameters of our study are modi- 
fied, we have tabulated the retention fraction for a rather 
comprehensive set of parameters associated with the selec- 
tion of primordial binaries and the behavior of mass trans- 
fer (Table 3). For Table 3, we have fixed the escape speed 
at fesc = 50 km s~^ and we have utilized the Maxwellian 
kick distribution with cr = 200 km s~^. The largest net 
retention fraction of ~ 8.3% (model 9) is only a factor 
~ 1.5 times larger than the retention fraction computed 
for the standard model (model 5). Thus, we can be secure 
that, for reasonable variations in the parameters listed in 
Table 3, the retention fraction is never much larger than 
the value computed for the standard model. 
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Table 3 

Retention Fractions for Various Parameter Sets 



Model 
Number 












n ■ f 
ycrlt 


Percent 
Retained 


1 


2.0 


-0.5 


1.5 


0.75 


1.0 


0.50 


3.23 


2 


2.0 


0.0 


1.5 


0.75 


1.0 


0.50 


5.39 


3 


2.0 


1.0 


1.5 


0.75 


1.0 


0.50 


7.87 


4 


2.5 


-0.5 


1.5 


0.75 


1.0 


0.50 


3.28 


5* 


2.5 


0.0 


1.5 


0.75 


1.0 


0.50 


5.79 


6 


2.5 


1.0 


1.5 


0.75 


1.0 


0.50 


8.03 


7 


3.0 


-0.5 


1.5 


0.75 


1.0 


0.50 


3.41 


8 


3.0 


0.0 


1.5 


0.75 


1.0 


0.50 


5.71 


9 


3.0 


1.0 


1.5 


0.75 


1.0 


0.50 


8.27 


10 .... 


2.5 


0.0 


1.0 


0.25 


0.3 


0.50 


3.70 


11 .... 


2.5 


0.0 


1.0 


0.25 


1.0 


0.50 


4.12 


12 .... 


2.5 


0.0 


1.0 


0.75 


0.3 


0.50 


4.20 


13 .... 


2.5 


0.0 


1.0 


0.75 


1.0 


0.50 


3.70 


14 .... 


2.5 


0.0 


1.0 


1.00 


0.3 


0.50 


2.94 


15 .... 


2.5 


0.0 


1.0 


1.00 


1.0 


0.50 


3.32 


16 .... 


2.5 


0.0 


1.5 


0.25 


0.3 


0.50 


2.23 


17 .... 


2.5 


0.0 


1.5 


0.25 


1.0 


0.50 


2.63 


18 .... 


2.5 


0.0 


1.5 


0.75 


0.3 


0.50 


5.15 


19 .... 


2.5 


0.0 


1.5 


1.00 


0.3 


0.50 


2.94 


20 .... 


2.5 


0.0 


1.5 


1.00 


1.0 


0.50 


3.32 


21 .... 


2.5 


0.0 


2.0 


0.25 


0.3 


0.50 


0.72 


22 .... 


2.5 


0.0 


2.0 


0.25 


1.0 


0.50 


1.15 


23 .... 


2.5 


0.0 


2.0 


0.75 


0.3 


0.50 


5.54 


24 .... 


2.5 


0.0 


2.0 


0.75 


1.0 


0.50 


5.93 


25 .... 


2.5 


0.0 


2.0 


1.00 


0.3 


0.50 


2.94 


26 .... 


2.5 


0.0 


2.0 


1.00 


1.0 


0.50 


3.32 


27 .... 


2.5 


0.0 


1.5 


0.75 


1.0 


0.25 


7.06 


28 .... 


2.5 


0.0 


1.5 


0.75 


1.0 


0.75 


3.38 



'Standard model. 
"IMF exponent. 

''Exponent for mass ratio distribution. 
"^Angular momentum-loss parameter. 
''Mass capture fraction. 
"Common-envelope efficiency. 

'Critical mass ratio that separates stable and unstable mass 
transfer when the secondary is radiative. 

6.3. Modified Lorentzian Kicks 

We now repeat the exercise of the last section, but 
with the modified Lorentzian kick distribution proposed by 
Paczynski (1990). The quahtative feature of Paczyhski's 
distribution that distinguishes it from the Maxwehian is its 
finite vahie at vanishing kick speeds. Physically speaking, 
this feature is unrealistic, since various stochastic processes 
associated with core-collapse and the SN explosion are cer- 
tain to deliver some net impulse to the newly formed NS; 
the characteristic minimum kick speed is likely to be of 
order 10 km s^^ rather than, say, 10 m s~^. However, the 
possibility that some finite fraction of NSs receive kicks 
between 10 and 50 km s~^ is not necessarily unrealistic. 

Table 4 summarizes the results of applying the mod- 
ified Lorentzian distribution, for Vesc = 50 km s~^ and 
uo(km s^^) = {100, 300, 500}. Of particular significance is 
the large fraction of NSs born in case D binaries that are 
retained in the cluster 6% for vo = 500 km s~^). This 
simply indicates that the Paczynski distribution, even with 
a large mean speed, allows a high percentage of the NSs 
born in isolation to be retained in a cluster. In fact, for the 
values of vq used to generate Table 4, the retention fraction 
of single NSs dominates over the contribution from any of 



the other five cases of mass transfer. Figure 14 illustrates 
this point clearly; the percentage of NSs retained via bi- 
nary channels is never more than a factor of two over the 
retained percentage of NSs born from isolated progenitors. 
6.4. Spatial Distribution and the Cluster Potential 

Up to this point we have discussed those calculations 
where the nominal central escape speed has been applied 
to all stars and binaries in question. The combination of 
competitive gas accretion processes, stellar collisions, and 
dynamical mass segregation in the early phases of clus- 
ter development and star formation may lead to a cen- 
trally concentrated population of massive stars (see Bon- 
nell. Bate, & Zinnecker 1998). However, the spatial dis- 
tribution would certainly have been finite and the same 
escape speed would not have applied to all stars. We in- 
vestigate the possibility that massive stars and binaries are 
born within a finite spherical volume, with a gravitational 
potential that is appropriate for a young globular cluster. 

For simplicity, we suppose that all massive single stars 
and binaries are distributed uniformly within a spherical 
volume of radius R about the center of the cluster. Thus, 
the probability that an object is located within a spherical 
shell of radius r and thickness dr is simply 

p{r)dr = ^^. (18) 

Furthermore, at the time of the SN explosion, we assume 
that the single star or binary is at rest. This is reasonable, 
since it is expected that the characteristic speed of these 
massive objects is < 5 km s~^. 

The background of lower-mass objects in the cluster pro- 
vides the net gravitational potential well in which the mas- 
sive stars reside. This assumes that all the cluster stars 
were formed at essentially the same time or that the low- 
mass stars formed first. We adopt a Plummer model for 
the potential (e.g., Binney & Tremaine 1987), given by 

$ = (19) 

where is the total mass in background stars, r is the 
distance from the cluster center, and b is the "core" radius 
of the model. In dynamical models of globular cluster 
evolution that include of the effects of tidal mass loss, it 
is often assumed that any star that crosses a sphere of 
radius rt (the tidal radius) is lost from the cluster (e.g., 
Joshi, Nave, & Rasio 2001, and references therein). The 
tidal radius is a function of position in the Galaxy; at 
a few kiloparsecs from the Galactic center, rt is of order 
100 pc for a 10^ Mq cluster. The escape speed, Vescir), at 
a radius r is obtained from the energy relation 

^ ,(r) = $(rt) - $(r) . (20) 



2 esc V 



However, for the purposes of this investigation we assume 
that ri is sufficiently large in comparison to any relevant 
radius in the cluster that we may drop ^(r-t), so that 
i'ggj.(r) = — 2$(r). In this case, the core radius, b, is a 
simple function of the central escape speed, Wesc(O): 

, 2GM, 



(21) 



The output of our population synthesis code is a two- 
dimensional grid of retention fractions as a function of the 
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Table 4 

Branching Percentages for Modified Lorentzian Kicks 



Outcome A B, C, D 



Total 5.43 21.44 2.96 15.30 6.76 48.11 

SlabloMT* 2.70 10.78 0.00 7.76 0.00 0.00 

Unstable MT (CE"^) 2.73 10.66 2.96 7.54 6.76 0.00 

Merger Following MT 3.22 10.64 1.94 5.86 0.33 0.00 



vo = 500 km s-i 

Unbound following SN 0.67 5.78 0.38 6.06 3.12 47.09 

Bound following SN 1.53 5.02 0.65 3.39 3.31 1.02 

Merger following SN 0.33 0.52 0.16 0.30 0.73 0.00 

Retained single NS 0.00 0.12 0.00 0.21 0.05 5.23 

Retained binary NS 1.06 4.29 0.44 2.51 1.81 1.02 

Total retained 1.06 4.40 0.45 2.71 1.86 6.25 



Wo = 300 km a-^ 

Unbound following SN 0.36 4.17 0.22 4.73 2.10 46.46 

Bound following SN 1.84 6.63 0.81 4.71 4.33 1.66 

Merger following SN 0.25 0.48 0.16 0.25 0.71 0.00 

Retained single NS 0.01 0.19 0.00 0.24 0.05 8.61 

Retained binary NS 1.36 5.90 0.63 3.69 2.66 1.66 

Total retained 1.37 6.09 0.63 3.93 2.71 10.26 



Vo = 100 km 

Unbound following SN 0.04 1.24 0.03 1.91 0.52 43.35 

Bound following SN 2.16 9.57 1.00 7.53 5.91 4.76 

Merger following SN 0.06 0.14 0.05 0.08 0.47 0.00 

Retained single NS 0.00 0.19 0.00 0.43 0.04 22.00 

Retained binary NS 1.66 8.87 0.93 6.35 4.52 4.76 

Total retained 1.67 9.07 0.93 6.78 4.56 26.76 



=MT = Mass Transfer. 
''CE = Common-Envelope. 



escape speed and the kick speed (see § 6.1). Combining 

eqs. (18), (19), and (20; dropping $(rt)), we obtain a dis- 
tribution in escape speeds for the spherical distribution of 
massive stars, 

p(Uesc) dUesc = 6 (^^^ Ugsc (1 " "^.J^/^ rfUesc , (22) 

where Wcsc = i^csc/i^osc(0) is a dimensionless escape speed. 

It is a simple matter to convolve the grid of retention 
fractions with both the distribution of kick speeds and the 
distribution of escape speeds to obtain a net retention frac- 
tion, as a function of Uosc(O) and R. We have tabulated 
(Table 5) the net retention fraction for different values of 
■1^680(0) and R using the MaxwcUian kick distribution with 
(T = 200 km s~^ and a cluster mass of = 10^. For 
t;esc(0) = 50 km s~^ and i? = 20 pc, the percentage of NS 
retained in the cluster is reduced by a factor of ~ 3 below 
the standard model value of ^ 5.6% (see Table 3). Thus, 
a realistic cluster potential and finite spatial distribution 
of stars is indeed an important consideration. 

6.5. Binary Evolution After the First Supernova 

The standard model (§ 6.2.1) has a very striking feature: 
massive secondaries (M2 > 10 Mq) are prevalent among 
the retained binaries following the first SN (see Figs. 11 
and 12). The majority of these massive systems have pe- 
riastron separations < 1 AU, which implies that most of 
the secondaries will begin to transfer material to the NS 
at some point. Furthermore, the extreme mass ratios sug- 
gest that the mass transfer will be dynamically unstable. 



resulting in a spiral-in of the NS into the envelope of the 

secondary. It should be noted that the evolution of the sec- 
ondary following mass transfer may not precisely resemble 
the evolution of an isolated star of the same mass (e.g., 
Braun & Langer 1995; Wellstein, Langer, & Braun 2001); 
in fact, the evolution may be qualitatively different. 

Before we discuss the possible outcomes of the spiral-in, 
we mention an important caveat. Extreme accretion rates 
(> 10~^ Mq yr~^) onto the NS - rates far exceeding the 
standard, radiative Eddington limit of ^ 10^^ Mq yr~^ - 
may be possible if the gravitational energy is lost to neu- 
trinos (e.g.. Chevalier 1993, 1996; Fryer, Benz, & Herant 
1996; Brown, Lee, & Bethe 2000). If this "hypercritical 
accretion" occurs while the NS spirals into the envelope of 
a massive secondary, it is likely that the NS will collapse 
into a black hole, although the three-dimensional nature 
of the hydrodynainical problem implies that this process 
is very uncertain. Obviously, this outcome is not desirable 
in regard to the retention problem, since a NS is lost if it 
is transformed into a black hole. We now proceed under 
the assumption that the NS does not undergo hypercriti- 
cal accretion during the spiral-in phase; however, the NS 
may still collapse to a black hole at a later stage. 

The envelope of the massive secondary may be success- 
fully ejected if the circularized orbital separation is > 1 AU 
(see Taam, Bodcnheimcr, & Ostriker 1978). This applies 
to only a few percent of the systems in Figs. 11 and 12. 
If the envelope is ejected, and the helium core of the sec- 
ondary is exposed, the formation of a second NS is possi- 
ble. However, the eventual SN explosion of the helium star 
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Table 5 

Modified Retention Fractions for a Finite Spatial Distribution of Stars 



«esc(0) 


R 


b 


Percent 


Reduction 


(km s-1) 


(pc) 


(pc) 


Retained 


Factor* 


30 


5 


9.88 


2.44 


0.94 


30 


10 


9.88 


2.10 


0.81 


30 


20 


9.88 


1.47 


0.56 


50 


5 


3.56 


4.44 


0.79 


50 


10 


3.56 


3.12 


0.55 


50 


20 


3.56 


1.81 


0.32 


70 


5 


1.81 


5.19 


0.64 


70 


10 


1.81 


3.38 


0.42 


70 


20 


1.81 


1.87 


0.23 



'Factor by which retention fraction is reduced 
below value obtained with R = pc. 



is likely to send the first- and second-formed NSs speed- 
ing out of the cluster, even if the kick to the second NS 
is small and the binary remains bound after the explosion 
(see § 5). 

The much more likely outcome among the retained NS 
binaries is a complete coalescence of the NS and the mas- 
sive secondary, resulting in the formation of a Thorne- 
Zytkow object (TZO; Thorne & Zytkow 1975, 1977; Biehle 
1991; Cannon 1993; Podsiadlowski, Cannon, & Rees 1995), 
where hydrostatic support is provided by gravitational en- 
ergy release or exotic nuclear burning processes near the 
surface of the NS. The ultimate fate of the NS is unclear. 
If the NS survives, it will probably emerge as a rapidly ro- 
tating object with the slow speed of the retained post-SN 
binary. However, it is possible, and perhaps likely, that the 
NS will collapse into a black hole during the late stage of 
massive TZO evolution (Podsiadlowski, Cannon, & Rees 
1995; Fryer, Benz, & Hcrant 1996). 

In addition to the high-mass systems in Figs. 11 and 
12. roughly 10% of the retained binaries have low- to 
intermediate-mass secondaries (M2 < 8 Mq), all of which 
are the product of dynamically unstable mass transfer 
prior to the SN. Mass transfer onto the NS in the cir- 
cularized binary is likely to be dynamically unstable for 
M2 > 4 M© (Podsiadlowski, Rappaport, & Pfahl 2001), 
while for secondaries of lower mass the system will exist 
for some time as a low- or intermediate-mass X-ray binary, 
which may ultimately yield a millisecond radio pulsar with 
a very low-mass companion. 

7. CONCLUSIONS 

The NS retention fraction calculated within our stan- 
dard model is ~ 5% for NSs born in binary systems. 

Reasonable variations of the parameters that describe the 
primordial binary population and binary evolution during 
mass transfer give a retained percentage between ^ 1 and 
8% (Table 3). If we distribute the massive binaries within 
a sphere of some finite radius and embed the population 
in a realistic background gravitational potential, the re- 
tention fraction may be reduced by a factor of ~ 2 — 3 
(Table 5). If we suppose that one- half of the massive stars 
in a young cluster are in binaries, then the retention frac- 
tion is further reduced by a factor of two. Therefore, a 
more realistic net retention fraction is probably of order 
1%, when we apply the Maxwellian kick distribution with 
a = 200 km s'^ 



As compared to the contribution from single stars, bi- 
nary systems do provide a much more efficient channel for 
retaining NSs when the characteristic kick speed is large. 
However, it appears the net NS retention fraction may not 
be sufficient to explain the abundance of NSs in globular 
clusters. In fact, even if as many as 10^ NSs are formed out 
of 10® stars, our standard model, combined with the binary 
fraction and realistic cluster potential, predicts that only 
100 NSs could have been retained. It is unlikely that such 
a small number of NSs is compatible with what is observed 
in certain clusters (e.g., 47 Tuc). We suggest that bina- 
ries alone to not provide a robust solution to the retention 
problem, and we now discuss alternative hypotheses. 



8. DISCUSSION 

Our standard model for the formation and retention 
of NS in globular clusters predicts a retention fraction of 

5%. Additional layers of realism, including a finite spa- 
tial distribution of stars embedded in a background clus- 
ter potential, and a binary fraction of massive stars, may 
reduce this number by more than a factor of four. We 
suggest that our standard model requires significant mod- 
ification in order for the results to be consistent with obser- 
vations. Here, we briefly speculate on possible alternative 
solutions to the retention problem. 

8.1. The Kick Distribution 

A very simple "fix" to the retention problem is to as- 
sume that the true underlying distribution in NS kicks 
has a much lower mean speed than the Galactic pulsars 
suggest. Such a distribution would require a substantial 
number of slowly-moving pulsars, which, for some reason, 
have not yet been detected. Maybe the pulsar sample 
is too small, or perhaps some systematic effect is unac- 
counted for in studies of pulsar kinematics, or both. A 
complete consistency check is difficult, if not impossible, 
and must incorporate all of the uncertain theory of single 
star and binary stellar evolution. In addition, we require 
some rudimentary understanding of the physical mecha- 
nism that produces the largest NS kicks. One possibility 
is that small kicks preferentially occur in binary systems, 
in which case the associated NS is likely to remain bound 
to its companion following the SN and thus will not appear 
as an isolated pulsar. 
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8.2. Implications of Wide, Nearly Circular High-Mass 
X-ray Binaries 

Recently, a new class of high-mass X-ray binaries 
(HMXBs) has emerged. These systems exhibit relatively 
low eccentricities (e < 0.2) and orbital periods sufficiently 
long (Porb ^ 100 d) that tidal circularization could not 
have played a significant role. A detailed analysis of one 
such system (X Per/4U 0352-^30; Delgado-Martf ct al. 
2001) revealed that the observed orbit is entirely consistent 
with the NS having been born with no kick whatsoever. 
Tabic 6 lists names and orbital parameters of systems that 
belong to the class of wide, nearly circular HMXBs. 

It is interesting to note that the present orbit of XTE 
J1543-568 (in 't Zand, Corbet, & Marshall 2001) is not 
precisely consistent with both zero kick and a perfectly 
circular prc-SN orbit; the eccentricity is too small by a 
factor of 2 — 3. One or the other assumption must be re- 
laxed. If we demand that the pre-SN orbit was circular, 
then it is straightforward to show that a very small kick 
of < 10 km s~^ is required to produce the observed eccen- 
tricity, and that the direction of the kick does not have to 
be finely tuned. 

With some perspective, we can motivate a phenomeno- 
logical picture that accounts for the population of wide, 
nearly circular HMXBs, and which is consistent with what 
we know about the Galactic NS population. Our model 
must satisfy three basic constraints. First, the systems 
listed in Table 6 are quite wide and thus have not ex- 
perienced the dramatic orbital shrinkage associated with 
dynamically unstable mass transfer; the mass transfer in 
these systems was likely in accord with the stable case Bg 
or Ce scenario described in § 4.2. We propose that a signif- 
icant fraction of those NSs whose progenitors underwent 
case Be or Ce mass transfer received only a small kick 
(e.g., < 50 km s~^). Second, the model should be able to 
approximately reproduce the numbers and properties of 
the observed population of bright low-mass X-ray binaries 
(LMXBs) in the Galaxy. Third, our basic picture should 
also be consistent with the observed kinematical distribu- 
tion of isolated pulsars in the Galaxy, on which the NS 
kick distributions are based. These latter two constraints 
are satisfied if we suppose that a NS received the usual 
large kick if its progenitor was allowed to evolve into a 
red supergiant (i.e., a single progenitor or case B/, C/, or 
D for a binary system). The standard formation channel 
for LMXBs involves a common-envelope phase in the case 
B; or C/ scenario (e.g., Kalogera & Webbink 1998), and 
so the NSs in these systems would have received the con- 
ventional large kicks by our hypothesis. Also, within this 
framework, isolated, fast-moving pulsars are likely to have 
come from single progenitors or wide binaries that were 
disrupted owing to the SN. 

We repeated the retention study with the following as- 
sumptions regarding NS kicks. If the orbit of the primor- 
dial binary is sufficiently wide that mass transfer begins 
when the star is a red supergiant (i.e., case B;, C;, or 
D), the NS kick is chosen from a Maxwellian distribu- 
tion with a = 200 km s~^, as in the standard model (see 
§ 6.2.1). However, if the mass transfer is stable and case 
Be or Ce, we utilize a Maxwellian kick distribution with 
(7 = 10 km s~^. Using the standard-model parameters 
given in Table 3, model 5, we calculate a net retention 



fraction contributed by binary channels of ^ 20%. There- 
fore, we see that the above scenario greatly reduces the 
severity of the retention problem. Figure 15 shows the 
distributions in orbital parameters and speeds for the re- 
tained NS binaries. On average, the systems in Fig. 15 are 
wider, more circular, and are moving more slowly than the 
systems in Fig. 11, which shows the same distributions, 
but for the standard model. 
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Fig. 15. — Distributions of binary parameters of systems tliat 
have undergone case B or C mass transfer, have been left bound 
following the supernova explosion, and have been retained in the 
cluster, where we have assumed that case Be and Ce receive kicks 
drawn from a Maxwellian with it = 10 km s~^. Compare this figure 
to Fig. 11. 

There is a plausible physical argument that may sup- 
port the empirically-motivated phenomenological hypoth- 
esis outlined above. Young, isolated, massive stars are 
observed to rotate at 20 — 50% of their breakup rates 
(e.g., Fukuda 1982). During the hydrogen-burning main 
sequence, the structure of the star changes sufficiently 
slowly that various hydro dynamical and magnetohydrody- 
namical processes should be effective in enforcing nearly 
uniform rotation throughout much of the star (Spruit & 
Phinney 1998; Heger, Langer, & Woosley 2000). Imme- 
diately following the depletion of hydrogen in the core, 
the structure of the star changes dramatically; the enve- 
lope expands to giant dimensions on a thermal timescale 
(- 10^ - lO'^ yr) and the core contracts while conserving 
angular momentum. If the helium core is exposed during 
this rapid expansion phase as a result of mass transfer in 
a binary system (case Be or possibly Ce; see Figs. 2 and 
3), the nascent helium star is likely to be rapidly rotating. 
On the other hand, if the star is allowed to evolve into a 
red supergiant. Maxwell stresses may strongly couple the 
mostly convective and very slowly rotating envelope to the 
core, causing the core to spin down to the angular velocity 
of the envelope (Spruit & Phinney 1998; Spruit 1998). 

This argument suggests a possible dichotomy in core- 
collapse dynamics between isolated stars and some stars 
born in close binary systems. It may be that a helium core 
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Table 6 

Orbital Parameters for Nearly Circular High-mass X-ray Binaries 



Object Porb(days) e /x(M)(Mq)^ Reference 



X Per/4U 0352+30 . . 249.90 ± 0.50 0.111 ± 0.018 1.61 ± 0.06 Delgado-Marti' et al. 2001 

7 Cas'' 203.59 ± 0.29 0.260 ± 0.035 ■ ■ ■ Harmanec et al. 2001 

XTE J1543-568 75.56 ± 0.25 < 0.03 8.2 ± 0.5 in 't Zand, Corbet, & Marshall 2001 

2S 1553-542'^ 30.60 ± 2.20 < 0.09 5.0 ± 2.1 Kelley, Rappaport, & Ayasli 1983 

GS 0834-430 105.80 ± 0.40 < 0.17 0.2 ± 0.3 Wilson et al. 1997 



"Mass function from X-ray timing. 

"^Orbital parameters determined from optical light curve. 

■^This orbital period is sufficiently short that tidal interactions may have circularized the orbit somewhat (Kelley, 
Rappaport, Ayasli 1983). 



exposed following case Be or Ce mass transfer is rotating 
much faster than the core of an isolated star at a late stage 
of its evolution. Dynamically, the collapse of a rapidly ro- 
tating core is certainly more phenomenologically complex 
than the collapse of an initially static core. However, it is 
not obvious a priori whether a rapidly rotating or slowly 
rotating core should ultimately yield a larger average na- 
tal kick to the NS. Perhaps the rotation axis provides a 
preferred direction for the escape of neutrinos or the for- 
mation of jets and the NS receives a kick perpendicular 
to the orbital plane. Another possibility is that rapid ro- 
tation stalls the collapse somewhat (e.g.. Fryer & Heger 
2000), allowing many rotations before the NS is formed. 
A large number of rotations of the collapsing core may have 
the effect of averaging out the asymmetries that give rise 
to large NS kicks (Spruit & Phinney 1998). These spec- 
ulations aside, we are motivated by empirical evidence to 
suggest that helium stars exposed following stable case Be 
or (possibly) Ce mass transfer produce NSs with small na- 
tal kicks, while NSs formed following mass transfer at a 
later stage of evolution may receive the conventional large 
kicks. 

8.3. Accretion- Induced Collapse 

Thus far, we have considered only massive stellar pro- 
genitors of NSs. However, if the mass of a white dwarf 
can be increased to the critical Chandrasekhar value (~ 
1.4 Mq), the white dwarf may collapse to form a NS. This 
accretion-induced collapse (AIC) scenario was proposed by 
Grindlay (1987; see also Bailyn & Grindlay 1990) in the 
context of globular clusters to explain a number of things 
regarding cluster NS populations, the retention problem 
among them. Two fundamentally different scenarios have 
been proposed for the formation of a NS via AIC, which 
we now discuss in turn. 

A white dwarf may be "grown" to the Chandrasekhar 
mass by the slow accumulation of material accreted from 
a stellar companion. In this scenario, if the white dwarf 
has a C/0 composition it is more likely to explode in a 
Type la SN, than to collapse to form a neutron star (e.g., 
Nomoto 1987; Rappaport, DiStefano, & Smith 1994). 
More favorable candidates for AIC are white dwarfs with 
an 0/Ne/Mg composition (e.g., Nomoto 1987; Nomoto & 
Kondo 1991). These are relatively rare white dwarfs with 
masses above ~ 1.2 Af©. To grow the white dwarf to the 
Chandrasekhar mass, however, rc;quires some fine tuning 
in the accretion rate, which must lie in the relatively nar- 
row range of about 3 - 7 x 10"'^ Mq yr"^ (Iben 1982; 



Nomoto 1982; Rappaport, DiStefano, & Smith 1994). For 

significantly lower accretion rates, the burning of hydrogen 
to helium is likely to be unstable, leading to hydrodynam- 
ical nova explosions, which may eject at least as much 
mass as was accreted (e.g., Prialnik & Kovetz 1995). For 
larger transfer rates the white dwarf atmosphere will tend 
to swell to giant dimensions and may overflow its Roche 
lobe, thereby losing much of the accreted matter. Per- 
haps the most promising cases for obtaining mass trans- 
fer rates within the above narrow range occur for thermal 
timescale mass transfer via the Roche lobe overflow in bi- 
naries with relatively unevolved companions in the mass 
range 1.5 — 2.5M0 (e.g., supersoft X-ray sources; van den 
Heuvel et al. 1992; Rappaport, DiStefano, & Smith 1994). 
Another possibility occurs for the case of accretion from 
the strong stellar wind of a low-mass giant (Iben & Tu- 
tukov 1984; Hachisu, Kato, & Nomoto 1999). 

A second possible AIC channel involves the coalescence 
of two white dwarfs (e.g., Nomoto 1987; Chen & Leonard 
1993; Rasio & Shapiro 1995). Two white dwarfs in a close 
binary system will be drawn together as gravitational ra- 
diation removes orbital angular momentum. The negative 
mass-radius exponent of a white dwarf implies that the 
less massive component will first fill its Roche lobe. Like 
the standard AIC scenario, which involves only one white 
dwarf, the double white dwarf merger model has been pro- 
posed as an evolutionary pathway to the formation of Type 
la SNe (e.g., Iben & Tutukov 1984; Webbink 1984; Saffer, 
Livio, & Yungelson 1998, and references therein); how- 
ever, it is now considered more likely that this will lead to 
disruption of the lighter white dwarf (see Nomoto & Iben 
1985). Under the assumption that the merger does not 
produce a Type la SN, in order for the double white dwarf 
system to ultimately yield a NS, the sum of the masses 
must exceed the Chandrasekhar mass. Approximately 1 
out of every 1000 primordial binaries should produce such 
a massive double white dwarf close enough to merge within 
a Hubble time (e.g., Han 1998; Nelemans et al. 2001). If 
a sizable fraction of these systems can collapse to form a 
NS, rather than explode as a Type la SN, then perhaps 
as many as 1000 NSs can be formed in this way in a glob- 
ular cluster. However, binary population synthesis in a 
globular cluster is substantially more complex than in the 
Galactic plane, owing to the dynamical interactions among 
binaries and single stars. 

We have not attempted to follow any of these channels 
in this work. This would certainly be a worthwhile future 
study to help quantify the formation rates of NSs via AIC 
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in both globular clusters and in the Galaetie plane. Fi- 
nally, we note that there is no obvious reason why a NS 
formed via AIC would be less subject to the same type of 
accelerations as those formed from collapsed cores of mas- 
sive stars, especially if the NS velocities arc acquired as a 
result of asymmetric neutrino emission or the post-natal 
electromagnetic rocket mechanism (Harrison & Tademaru 
1975; see also Lai, Chernoff, & Cordes 2001). 

8.4. Supermassive Globular Clusters 

Tidal stripping of globular clusters is a theoretically 
well-studied phenomenon (e.g., Chernoff fc Weinberg 1990; 
Takahashi & Portegies Zwart 2000; Joshi, Nave, & Rasio 
2001). It has been shown (see Joshi, Nave, & Rasio for a 
recent discussion) that, for a range in parameters that de- 
scribe the initial cluster equilibrium model, a cluster may 
disrupt in the tidal field of the Galaxy in less than 10^0 yr, 
owing to the combined effects of mass loss during stellar 
evolution and the diffusion of stars across the cluster's tidal 
boundary (effectively, its Roche lobe). The survivability 
of a cluster depends on its location in the Galaxy, its cen- 
tral concentration, and on the shape of the cluster IMF (in 
particular, the slope of the IMF above ^ 2 Mq). Clusters 
with a high central concentration and a small proportion 
of massive stars are more likely to survive to core collapse, 
although as much as 90% of the initial mass of the cluster 
may be lost before this phase is reached (Joshi, Nave, & 
Rasio 2001). 

The idea that a cluster may lose a very significant frac- 
tion of its mass, but still "survive" in some sense, brings 
to light the interesting and very real possibility that some 



of the globular clusters that presently have a mass of 
~ 10^ Mq are, in fact, remnants of clusters with an initial 
mass in stars of > 10^ Mq. At least one such supermas- 
sive cluster has been discovered in the Andromeda galaxy 
(Meylan et al. 2001), and it has been speculated that this 
cluster is actually the core of a dwarf elliptical galaxy. For 
a cluster of initial mass 10^ Mq, a net NS retention frac- 
tion of a few percent, along with a standard IMF, implies 
possibly a thousand NSs at the current epoch, which may 
be sufficient to explain the present pulsar and X-ray binary 
populations in globular clusters. 

We are not the first to consider this rather extreme 
possibility. Motivated by the hundred-fold overabundance 
(per unit mass) of bright X-ray sources in globular clus- 
ters relative to the Galactic disk, Katz (1983) suggested 
that perhaps some clusters lose all but ~ 1% of their mass 
through evaporative processes. Although it is now be- 
lieved that 3- and 4-body dynamical scenarios (e.g., Rasio, 
Pfahl, & Rappaport 2000), and perhaps tidal capture (e.g., 
Fabian, Pringle, & Rees 1975; DiStefano & Rappaport 
1992; Podsiadlowski, Rappaport, & Pfahl 2001), can ex- 
plain the overabundance of X-ray binaries, excessive mass 
loss from initially supermassive globular clusters is still an 
interesting possibility for explaining the large numbers of 
NSs found in clusters today. 
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APPENDIX 

A. ESTIMATE OF THE BINARY RETENTION FRACTION 

A semi-analytic approach is used to compute the probability that a binary is bound following the SN explosion and is 
retained in a cluster with a given escape speed. For an appropriate choice of pre-SN orbital separation and component 
masses, the calculated retention probability is a fair estimate of the net NS retention fraction contributed by all binary 
stellar evolution channels (see § 5). 

In what follows, we consider a circular pre-SN binary of total mass Mt and with separation o, so that the relative 
orbital speed is Worb = [GMb/ a)^^"^ . We assume that the explosion is instantaneous and leaves a neutron star remnant of 
mass Mns- Furthermore, we neglect the effect of the SN ejecta on the secondary. The kick speed imparted to the NS is 
Vk and the systemic mass after the explosion is M^. 

It is useful here to introduce a set of dimensionless variables. All speeds are expressed in units of the pre-SN relative 
orbital speed, ?;orb, and are denoted by the variable w with an appropriate subscript (e.g., Wk = Wfe/'^orb is the dimensionless 
kick speed). The fractional mass loss in the explosion is given by A = 1 — Af^/Aft. In place of the secondary mass M2, we 
use the post-SN mass ratio q' = Af2/AfNS- Finally, the variable u represents the cosine of the angle between the direction 
of the kick and the direction of the pre-SN relative orbital velocity; note u = when the kick is directed perpendicular to 
the orbital plane. 

The orbital energy, E' , of the post-SN binary is proportional to the dimensionless quantity (e.g., Hills 1983; Brandt & 
Podsiadlowski 1995) 

E' -1 + 2/:^ + wl+2uwk . (Al) 
For the purposes of § 5 it is sufficient to consider A < 1/2, in which case there is ci mininiuni kick, W]^ min ? 

required to 

unbind the system, realized when u = 1: 

W^fc,min = -l+[2(l-A)]l/2 . (A2) 

Likewise, for there to be a finite probability that the system remains bound, Wk must be less than some large value, 
■w^fe.max, corresponding to E' = Q and u = — 1 in eq. (Al): 

U^fe,ma. = l+[2(l-A)]l/2 . (A3) 
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If Wk > Wk,max, the system is guaranteed to be disrupted. For a given Wk,min < Wk < Wk^max and A < 1/2, u must be 
less than a maximum value, Umax, for the binary to remain bound following the explosion: 

u<u^^^ = -^{l-2A-wl) . (A4) 

2wk 

When the kick speed is large {wk > 1), we see that timax < 0. Therefore, if the directions of the kicks are preferentially 
aligned perpendicularly to the orbital plane (i.e., u ^ 0), rather than distributed isotropically, wc expect somewhat fewer 
bound and retained systems for a mean kick speed that is larger than the typical orbital speed (see Brandt & Podsiadlowski 
1995). Conversely, if the mean kick speed is small, perpendicular kicks tend to yield an increase in the number of retained 
binaries, although the baseline retention fraction (for isotropic kicks) is also larger in this case. 

The center-of-mass (CM) speed of a bound binary determines whether or not the system will be retained in the cluster. 
The CM speed, wcm, following the explosion is given by 

wcM = ^[{q'Af - 2q'AuWk + wlY/^ . (A5) 

If Wesc is the dimensionless central escape speed of the cluster, then the probability that a bound post-SN binary is 
retained is simply a step fimction, S{wcsc — 'U'cm), equal to imity for Wcsc — wqm > and vanishing otherwise. Taking A, 
q' , and w^sc to be fixed parameters, we obtain the retention probability, Pr, as a function of Wk by integrating over u: 

/min(l,Mniax) 
dup{u) S{wesc - Wcm) , (A6) 

where min(l,Minax) is the minimum of 1 and Umax- For isotropically distributed kick directions, the distribution fimction 
for u is simply p{u) = 1/2. Convolution of Pr with the distribution in dimensionless kick speeds, p{wk), yields the total 
probability, Pr.tot, that a bound binary is retained in the cluster after the SN: 

-Pr,tot(A,g',U;esc) = / dwk p{wk) Priwk] A, q' , Wesc) (A7) 

Jo 

B. SUPERNOVAE IN ECCENTRIC BINARIES 

In this Appendix, wc present a flexible, computationally convenient formulation of the equations that describe a binary 
system following an asymmetric SN (SN) explosion of one of the components. We allow for the possibility that the pre-SN 
binary is eccentric, and we consider the effects of instantaneous mass loss from the exploding star and an impulsive kick 
delivered to the newly- formed compact remnant. Also included in our analysis is the effect of the SN blast wave on the 
companion to the exploding star. Furthermore, if the binary is disrupted following the SN, we calculate the asymptotic 
velocities of the components. Our approach differs from previous studies (e.g.. Hills 1983; Brandt & Podsiadlowski 1995; 
Tauris & Takens 1998) in that we use mathematically compact vector expressions to describe the binary system after the 
explosion. It is straightforward to directly implement this vector formalism in a computer code, since vector arithmetic 
can be performed using simple array operations. 

Consider a pre-SN binary system that consists of stars with masses Mi and M2 in an orbit with semimajor axis a and 
eccentricity e. The Keplerian orbital frequency is given by n = (GMb/a^)!/^, where Mb = Mi + M2. Relative to the 
center-of-mass (CM), the positions of the two stars at some time t are ri{t) and r2{t), and the corresponding velocities 
are vi{t) and V2{t). The relative positions and velocities are given by r{t) = ri{t) — r2{t) and v{t) = vi{t) — V2{t), 
respectively. 

It is convenient to have a coordinate-independent description of the binary system. Such a description is provided 
by the two conserved vectors of the Kepler problem, namely the angular momentum per unit reduced mass, h, and the 
Laplace-Runge-Lenz (LRL) vector, e (e.g., Goldstein 1980; Eggleton 1999): 

I. vxh r 

Note that h points perpendicular to the orbital plane and has a magnitude h = fla^^l — e^)^^"^, and e points in the direction 
of periastron of star 1 and has a magnitude equal to the orbital eccentricity, e. By convention, boldfaced characters denote 
vectors, while the same characters with normal typeface denote the magnitudes of those vectors. 

Since we allow for the possibility that the pre-SN binary is eccentric, we must take some care in computing r and v at 
the time of the explosion. We assume that there is no preferred position along the orbit for the explosion to take place; 
therefore, the explosion probability per unit time is constant. No closed form expressions exist for r and v as functions of 
time, and so we must be content with a parametric representation. Consider a Cartesian coordinate system with x-, y-, 
and 2;-axes defined by the directions of e, h x e, and h, respectively. In terms of the eccentric anomaly, E, the dynamical 
equations read 

ntp^E-esinE (B2) 
a; = a(cos£;-e) ; y = a{l - e^)^/^ sinE (B3) 

v^ = -^smE ; Uy = ^(1 - e^)^/^ cosiJ , (B4) 
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where tp is the time elapsed since periastron passage. With a randomly selected value for flip, the corresponding value of 
E is obtained by solving eq. (B2) numerically, and the relative position and velocity vectors can be readily computed. 

At the randomly selected time, star 1 undergoes a SN explosion. We assume that the explosion is an impulsive event, 
meaning that the direct dynamical influence of the explosion occurs over a time that is much shorter than the orbital 
period. In other words, it is assumed the SN explosion and the associated blast wave have an instantaneous effect on the 
masses and velocities of the binary components. The envelope of star 1 is ejected, exposing a remnant of mass M[ with a 
new velocity v[ = vi + Avi, where Avi is the kick velocity. The magnitude and direction of the kick velocity are chosen 
from appropriate distributions (see § 3). After a negligibly short time (the time it takes the blast wave to cross the orbit), 
a small fraction of the blast wave will interact with star 2, resulting in a new mass and velocity v'2 = V2 + Av2, 
where Av2 is directed antiparallcl to r (see Wheeler, Lecar, & McKee 1975; Fryxell & Arnett 1981). If star 2 is still on 
the main sequence at the time of the explosion, it is expected the SN ejecta has only a small effect on star 2 and on the 
binary orbit. However, if star 2 is a giant at the time of the SN, a large fraction of its envelope may be stripped by the 
blast wave (see, however, Livne, Tuchman, & Wheeler 1992; Marietta, Burrows, & Fryxell 2000), and so we consider this 
possibility in our mathematical formalism. 

The combined efii'ects of mass loss and the velocity perturbations received by the binary components yield a CM velocity 

v'cM = ^[Miv[+M^v'2] 

M2 AMi Ml AM2 \ Ml Mi , , 



Ml Mb Ml Mb J Ml M^ 

where AMi = Mi — M[ is the mass of the ejected envelope of star 1, AM2 = M2 — M2 is the mass stripped and ablated 
from star 2 (Wheeler, Lecar, & McKee 1975), and Ml = M{ + M2. The orbital parameters following the explosion may 
differ dramatically from their initial values. In fact, the explosion may disrupt the binary entirely. The post-SN orbital 
parameters are determined by the new specific angular momentum, h' , and the new LRL vector, e': 

V = .x.' ; e' = l^-I^ (B6, 

The binary is gravitationally bound following the SN if e' < 1. In this case, the post-SN semimajor axis is given by 

a' = ^ . (B7) 

GMlil - e'2) ^ ^ 

It is sometimes interesting to know the spin-orbit misalignment angle, 7, of the compact remnant or its companion 
following the SN (e.g., Brandt & Podsiadlowski 1995; Kalogera 2000). Star 2 will have the same rotation sense as the 
orbit following a phase of mass accretion. In this case, the spin of star 2 preserves the direction of the pre-SN orbital 
angular momentum, and the cosine of the misalignment angle is simply 

cos J = h ■ h , (B8) 

where hats denote unit vectors. Star 1 may likewise spin in the direction of the orbit owing to tidal coupling; however, 
this is not necessarily true for the remnant (see Spruit fc Phinney 1998). 

On the other hand, if e' > 1, the compact remnant and star 2 are not gravitationally bound, and we would like to 
compute the asymptotic speeds of the components relative to the pre-SN CM velocity. In the new CM frame, the two 
objects recede along hyperbolic trajectories. As a function of the true anomaly (also the polar angle in the new orbital 
plane), 9, the relative separation increases according to r{0) oc (1 + e'cos^)"^. Clearly, r approaches infinity as cos^ 
approaches the value — 1/e'. For large r, the direction of the relative velocity is nearly radial, and so the relative velocity 
at infinity, is given by 



1 / 1 \ / 
-— e + i 1 2 I h xe 



(B9) 



where 



,^ = ^(e'^_ 1)1/2. (BIO) 

Given Voo and v'^^^, the asymptotic velocities of the components relative to the pre-SN CM velocity can be computed: 

M2 , M{ , , 

Vl,oc = J^'"oo + VcM ; '"2,00 = -Jf'"'^ + '"CM ■ (Bll) 
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